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ABSTRACT 



The effects of point defect implantation in copper and tungsten 
crystal lattice have been studied by computer simulation techniques. 
Vacancies, interstitials, and replacement impurities have been 
created in the first five layers of the free (100) surface of these 
crystals. The subsequent binding energies of these defects in 
tungsten were compared with experimental temperature dependent de- 
sorbtion peaks, corresponding to binding energies of neon defects 
in a tungsten crystal. Interstitial and replacement impurity 
positions in the first three to five layers were found that seem 
to correspond to the experimental data. Significant results were 
also obtained which were associated with general surface effects, 
especially c:rowdion migration. 
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I. INTRODUCTION 



Extensive research has taken place in the last decade in the 
area of computer simulation of radiation damage in crystal lat- 
tices. Two major areas of simulation have been defined. "Dynamic 
simulation" suggests the firing of an atom or ion against a crystal 
and the observation of the resulting many-body collisions. Ex- 
amples include sputtering simulation [l,2,3l, in which atoms are 
ejected from the surface of an ion-bombarded crystal; and chan- 
neling simulation t4] , in which ions are fired down open channels 
in non-close-packed structures such as body-centered cubic and 
diamond lattices. "Static simulation", on the other hand, is con- 
cerned with equilibrium positions and energies for point defects 
in crystals [5,6,71. This latter area was the concern of this 
research. Specifically, this simulation attempted to correlate 
equilibrium potential energies of point defects with experimentally 
determined binding energies of point defects in Tungsten [8,91. 

A. HISTORICAL RA.CKGROUND 
1 . The Problem 

Historically, all crystal dynamics computer simulation has 
been based on the assumption that the complicated many-body problem 
can be reduced to many two-body problems. This assumption has re- 
peatedly been shown to be a valid one when employed incrementally. 
Incremental calculations are necessary since a complete solution 
in closed form is impossible. Small time increments At approximated 
the true time differential dt of the impossible closed form 
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solution. Specif ically 5 the desired type, size, and orientation 
of crystal lattice is stored in the computer, appropriate inter- 
atomic potentials are chosen, and all mutual forces between atoms 
are calculated, based on the analytic potential functions. Point 
defects are introduced, and each atom is then allowed to move in- 
crementally, based on these forces and Newtonian Mechanics C 5l • 
Through proper choice of time increment duration, damping of forces 
and velocities in each time increment, and sufficient repetition 
of the procedure, realistic results are obtained. 

2 . The Pioneers and Their Contributions 

Pioneering work in this field began in the early 1960’s 
at the Brookhaven National Laboratory. Gibson, Goland, Milgram 
.and Vineyard (GGMV) [ 5l published the results of extensive work 
in both static and dynamic simulation. In static simulaticr. they 
determined equilibrium positions for interstitials and associated 
potential energies of formation. In dynamic simulation they in- 
vestigated momentum propagation directions of energetic knock-on 
atoms (focusing), collision chains, and related topics. They used 
a central-difference method to obtain velocities and positions from 
calculated forces. All work was done with copper, and results were 
correlated with experimental data. Johnson and Brown (JB) [6] did 
extensive work in static simulation, again with copper. They es- 
tablished that only one stable position exists for a single face- 
centered cubic (FCC) self-interstitial: the ( lOO) split inter- 
stitial. Johnson [lO] later published further work in this area, 
with formation and activation energies for various point defects. 
Enginsoy, Vineyard, and Englert (EVE) [?] and Johnson [llj 
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repeated most of the earlier calculations in GGMV and Johnson 
for the body-centered cubic (BCC) case, based on O' iron. They, 
too, established the existence of only one stable interstitial 
position, a (lio) split interstitial. 

Girifalco and Weizer (GF) [ 12 ] calculated Morse Potential 

$.. = D[exp {-20i:(r..-r )}- 2 exp{ -0!(r . .-r }] 

parameters for various metals, based on experimental values for 
the energy of vaporization, the lattice constant, and compressi- 
bility. Resulting elastic constants and equations of state agreed 
satisfactorily with experiment. Girifalco and Weizer [ 13 ] later 
published results of using these Morse parameters in simulating 
vacancy relaxation dynamics. Anderman [l4l used GW’s technique 
of calculating Morse parameters, but instead of summing over an 
entire crystal, (GW calculated out to the 150th nearest neighbor) 
Anderman found parameters as a result of summing out to second, 
third, and fourth nearest neighbors, for use in short-range approxi- 
mations . 

Harrison Cl, 2, 3 ] has investigated sputtering phenomenon and 
other surface effects with a modified Brookhaven model, the most 
significant change being the use of an average force method C 15 ] 
instead of the central difference method in integrating the 
equations of motion. (See Appendix C. ) He has also calculated 
repulsive potentials of the Born-Mayer type j ) ) 

for many combinations of atoms and ions based on secondary elec- 
tron emission, and Hartree-Fock atomic electron distributions I I 6 ] . 
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3 . The Potential Function Problem 

The most difficult problem encountered by computer sim- 
ulation has been the proper choice of the potential function. No 
simple analytic expression, based on either theory or experimental 
data, has ever been found that completely describes crystal dy- 
namics Cl?!, although many analytic expressions are partially cor- 
rect. The problem has been three-fold: First, present analytic 
expressions have narrow regions of validity, i.e., some correctly 
describe atomic behavior at equilibrium distances, but fail at 
shorter or greater distances. Second, some analytic expressions 
are limited because they only apply to interactions between iden- 
tical atoms. Third, the assumed functions have spherical symmetry, 
and are technically limited to interactions between closed shell 
atoms or ions Cl^l Although our assumption of a. spherically sym- 
metric potential in crystals is only approximately correct, it is 
nevertheless a very good approximation for FCC structures, and a 
reasonably good approximation for BCC structures. It is grossly 
in error when applied to diamond structures. 

The atomic potential, with the familiar potential well, 
sharply repulsive wall and gently attractive tail, varies greatly 
between different pairs of atoms; Since theory can give only ap- 
proximate parameters for this complete potential function, experi- 
mental data have been used extensively in the formulation of po- 
tentials. Other avenues have been opened by computer simulation. 
Since potential well depths are typically on the order of a few 
eV, the characteristics of the well can be ignored in high energy 
dynamic simulation. Low energy dynamic simulation, and even static 
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simulation, have also been based upon this approximation with 
useful results. Historically this is how crystal simulation 
began. GGMV C 5 ] employed a purely repulsive potential of the 
Born-Mayer (BT4) type and applied external forces on all crystal 
boundaries to hold the crystal together. JB L6] used basically 
the same technique. For improved equilibrium studies a potential 
with a well was necessary so GW C 32,131 used a Morse potential in 
their simulation. The Morse function, however, fails at strongly 
repulsive distances. To satisfy the need for a more versatile 
potential, capable of handling both high energy and near-equili- 
brium dynamics, composite potentials were developed, which re- 
semble BM or Bohr 



(v^. = - — exp(A+Br^ )) 

ij 

functions at short separations, and Morse functions at equili- 
brium and greater separations. Specifically, EVE [?] combined 
a screened Coulomb or Bohr potential, a BM potential, and a Morse 
potential, in the higher repulsive, lower repulsive, and attractive 
regions, respectively, of the atomic potential. 

Johnson [lO,ll] in his later papers, used three cubic equations to 
approximate the true potential. Anderman C l4l and Harrison Cl, 2, 
3 j 4] have used the BM repulsive term together with a Morse well 
and attractive tail, smoothly fit together by a cubic equation in 
the region near their intersection. 

4 . The Point Defect Problem 

All early simulation was done with homogeneous systems: 

All atoms were exactly the same, limiting high energy dynamic 
simulation to bombardment by atoms identical to the lattice atoms, 
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and static simulation to consideration of only the vacancy and 
self-interstitial cases. This limitation was forced by the po- 
tential functions, because parameters for the Morse function were 
based on experimental data for homogeneous media [ 123 . The 
methods used could not yield parameters for different-atom pairs. 

BM parameters, however, are obtainable for different-atom pairs, 
by methods such as the Hartree-Fock method Cl6] mentioned pre- 
viously . 

In spite of these limitations, dynamic computer simulation 
of bombardment by foreign atoms, or static simulation of foreign 
interstitials, can be done by two alternate methods. First, 

Harrison [4l has neglected the attractive interactions and has 
done foreign particle dynamics using repulsion only. Alternately, 
Johnson Cll] has derived a cubic equation for a complete potential 
with a potential well, based on limited experimental data on car- 
bon defects in iron. However, experimental substantiation for a 
foreign-particle potential well is much more difficult than for 
an identical atom potential well. 

B. THE EXPERIMENT 

The experimental data which this simulation proposed to explain, 
were published by Kornelsen and Sinka (KS) L8] . They have bom- 

4- + + + 

barded a clean (100) tungsten surface with Ne , Ar , Kr , and Xe , 
in the energy range of 40 eV to 5 keV . The subsequent '^damaged” 
crystal was heated at a constant rate, and gas desorbtion rates 
were measured. Instead of a constant desorbtion of ions, various 
distinct peaks were found, categorized into two basic types: a 
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single large peak at iSOO^K, the same for all four ions; and four 
or five smaller peaks in the 400 ^K to l650^K range, which were 
not in the same position for all four ions. These latter peaks 
were postulated to correspond to binding energies of various point 
defects in the first few layers of the tungsten crystal. (See 
Figure 4. ) 

This simulation used Harrison^s assumption of a repulsive 
potential only, for interactions between a foreign point defect 
and other atoms in the lattice. When investigating neon defects 
in tungsten, all tungsten lattice interactions were based on com- 
posite Morse and repulsive Born-Meyer potentials, and all neon- 
tungsten (Ne-W) interactions were based on a purely repulsive BM 
potent ial . 
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.II. OBJECTIVE 



The long-range objective of this simulation was to correlate 
simulated and experimental binding energies of neon point defects 
in tungsten. Since the assumption that all Ne-W interactions are 
purely repulsive was not realistic, the degree to which subsequent 
simulation results are valid must be based on a known standard. If 
the simulated binding energies are not correct, a valid correction 
factor can be applied, if derivable from the known standard. One 
standard which proved to yield this information was the tungsten- 
tungsten (W-W) interaction. A tungsten point defect could be 
treated as an atom of the lattice, and given an interatomic po- 
tential identical to all other lattice atoms, the composite Morse 
and BM potential. This was Method 1, A tungsten defect could 
also be treated as foreign, and allowed to interact with other 
lattice atoms with a repulsive potential only (Method 2). If a 
specific tungsten point defect is treated by both of these methods, 
an empirical relationship between the repulsive potential assumpt- 
ion and the ^^true” potential for W-W interactions is obtained. 

The objectives of this research were fourfold: 

1. Demonstrate that the two methods of treating tungsten point 
defect in a tungsten latt ice yield basically the same physical 
results, and agree with published results [7,n] concerning 
split interstitial positions. 

2. Develop a general empirical relationship between the binding 
energies derived by the two methods. 
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3. Obtain values for binding energies of neon defects in all 
possible positions in a tungsten surface. Transform these 
values to more realistic ones using the empirical relationship 
derived in 2. 

4 . Compare these results with KS^s experimental data. 
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Ill . THE MODEL 



A . THE CRYSTAL 

The model used in this research is the Gay-Harrison C 19 ] model, 
with modifications by Levy [20] . Johnson C21], Effron [22], and 
Moore C 23 I . Abbreviations in brackets refer to computer program 
names for the variable in question. 

Both copper and tungsten crystals were simulated. Copper was 
simulated only to provide an interface between this research and 
published simulation results. Copper forms a face-centered cubic 
crystal with an experimentally determined lattice constant, (LC) 
or cube edge distance of 3«6l5S. The lattice unit (LU), defined 
^LC. is I. 8075 X; and the nearest neighbor distance, as in all 
FCC structures, is \T2LU. Tungsten forms a body-centered cubic 
crystal, with a LC of 3*l6^, a LU of 1.58^, smd a nearest neighbor 
distance, peculiar to all BCC structures, of \[3 lU. All distances 
in the program are measured in LU. The program could construct 
(100), (110), and (111) orientations of face-centered and body- 

centered cubic structures. The copper crystal size was 8x8x8, 
and contained 256 atoms for the (100) orientation. 

The major portion of the simulation was done on the (100) 
orientation of tungsten, corresponding to KS^s experimental work. 
This tungsten crystal size for Neon point defects was 10 x 10 x 10, 
and contained 250 atoms. Some W-W simulation was done on a 
14 X 14 X 14 crystal; the reasons are explained in RESULT S . The 
bottom two layers of the lattice were not allowed to move, al- 
though they had potential energy, and exerted force on all atoms 
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in the crystal. The other eight layers were completely free to 
move, and were included in the dynamic calculations of each time- 
step. 

The surface layer (Y = 0) and the second layer (Y = 1) were 
moved forward, simulating actual surface relaxation in the crystal. 
This relaxation was calculated by Moore C 24 ] using simulation 
techniques, and tested against previous results by Burton and 
Jura L 25 ]. Definitions and use of mobile layers, relaxation, etc., 
were analagous in the copper model, as were all other aspects of 
the model to be described in this chapter. 



B. THE POTENTIALS 

1 . The W-W Composite Potential 

The attractive potential used was the Morse potential, with 
tungsten parameters calculated by GVv Cl2]. The interaction energy 
0^^. , of a pair of particles i and j is: 

0 . . dC exp{ -2a(r..-r )} - 2 exp[ -0'(r..-r )}] (1) 

ij ^ ' ij o^ ^ ' ij o' ' ' 



where D [dCON] is the dissociation energy of the pair, r^CPE] is 
the equilibrium separation, r^^.LDISTj is the actual separation, 
and OL [alpha] is a constant. 

The repulsive potential is of the BM type, with Harrison^s 
Hartree-Fock parameters. The interaction energy , is 



V. . 



exp 



(A+Br , .) 

ij 



( 2 ) 



where B [exb] is always negative, A [exa] is always positive, and 
r.. [dISTJ is the actual separation. The constants [eXA] and 
[exb] are peculiar to the V7-W interaction. 
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The ranges of the W-W composite potentials were as follows 
the BM repulsive potential operated from 0 to 1-5^; and the Morse 
potential from 2 % CrOEb] to 5.38^ CroEC] . In LU, the dimension 
in which all calculations were done, these constants v/ere . 9494 , 
1 . 2658 , and 3*4000 LU. CrOEcJ was chosen to include interactions 
out to the fourth nearest neighbor (NN4) at \Tll LU = 3-31? LU but 
not NN5 interactions, at ^ 12 LU = 3*464 LU. Note, however, that 
slight displacements of NN5 ’s might allow their inclusion in po- 
tential and force calculations. The gap between 1 . 5^ and 2^ was 
filled with a cubic function, which matched to the other po- 
tentials and slopes at [roEA] and [roEbI . 

2 . Purely Repulsive Potentials 

For foreign point defect interactions, i.e., Ne-W, or 
W-W Method 2, a repulsive potential only was used. The potential 
was again a BM, with the constants labeled CpEXA] and CpEXb] . For 
the W-W, Method 2 interaction, LpEXA] = [eXA“] and CpEXB] = [exbI . 
Ranges for foreign point defect interactions, however, were dif- 
ferent, and the potential itself was modified at the cutoff point. 
Whereas the BM part of the composite potential extended out to 
about .95 LU, the modified BM potential used for foreign defects 
was allowed to extend to \f3 LU, corresponding to the NNl distance 
[roe]. Cutting the potential off at [rOe] left a step of about 
.05 eV for Ne-W (.2 eV for V^/-W) at the NNl equilibrium position. 
Since neither discontinuities nor repulsive potentials were de- 
sired at this equilibrium position, we ’^eroded” LIS'! the potential 
by subtracting V(CR0e1 ) , or about .05 eV from V(r^^) for 
r. . < [roe 3 . Calculated forces, based on these eroded potentials 
must be modified also, but for a different reason. It is possible 
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to conceive of a case where an atom is further away from the 
defect than [roe] at the beginning of a timestep, but closer than 
[roe] at the end. The force has essentially ^turned on'’ in the 
middle of the timestep. The modification gives such an atom a 
force which is less than the final force by approximately a fac- 
tor proportional to the ratio of distance traveled outside CrOE] 
to the total distance traveled during the timestep Cl5l. (See 
Appendix B. ) 

C. THE TIMESTEP 

Motion caused by these forces must be found by an approximate 
numerical method of time integration. As in previous work with 
this model, the average force method [l5] was used. In this 
method, all mutual forces were calculated in subroutine STEP. 
Based on these forces, new temporary velocities and positions 
were found. Forces were again calculated, based on the temporary 
positions. The final positions were then calculated from the 
average of these two force determinations. All velocities were 
then either zeroed or halved as a damping method. This consti- 
tuted one timestep. 

For this average force method to work properly, the AtCur] 
which approximates dt in the integration must be kept small . Too 
small a value for [dt], however, would result in excessive com- 
puter time. The choice of [dtJ was also complicated by the fact 
that Cdt] must be kept much smaller earlier in the program, when 
velocities, forces, and energies are large, but can be allowed to 
grow larger as the simulation approaches equilibrium. For this 
reason, at the end of each timestep, a new LDt] was calculated 



19 



for use in the next timestep. The parameter chosen to control 
Cot] was CDTl],the distance, measured in LU, which the most 
energetic atom was allowed to move before starting a new timestep. 
LdTIJ has varied between .001 and .02, depending on such con- 
ditions as original position of the point defect, relative masses 
of atoms, etc. In general, [ DTll must be kept very small when 
high velocities are expected, and can be increased when all motion 
is expected to be ^^sluggish”. In actual practice [dt] and [dTI] 
are related to both the? velocity of each particle and the force 
on each particle. To insure that no particle traveled more than 
CdTiI we ensured that [ Dt] was small enough so that neither the 
velocity of the most energetic atom nor the force on the most 
.stressed atom would result in motion greater than LdTi]. (See 
Appendix C. ) 

D. FOREIGN INTERSTITIALS 

1 . Unequal Mass Implications 

Many changes in the program were necessary when a foreign 
defect was included in the lattice. These changes were especially 
necessary when the defect was much lighter than the lattice atoms; 
i.e., neon in a tungsten lattice. First, in the average force 
calculations, a separate section had to be added for the calcu- 
lations for the primary, or ^*bullet’^, based on the bullet mass 
CbMAS] . Second, the potential energy between two unequal mass 
atoms was split in proportion to their reduced masses (^ee 
Appendix D) . Third, the section that determined the new timestep 
duration was originally based on the lattice atom mass. Since the 
light interstitial is usually the most energetic or most stressed 
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atom, very erratic behavior was observed until the timestep du- 
ration calculations were revised to handle tivo different masses 



(see Appendix C) . Finally, a significant mass difference between 
defect and lattice atom required a reduction in [ 0 x 1 ]. For the 
neon-tungsten simulation, a [dTiI of .5% was used, 

2 . Ionization State and Repulsive Potentials 

The major portion of the Ne-W work was done with the , 

assumption that tungsten was in a +6 state, and neon was neutral 
in the lattice. Experimentally KS fired neon in a +1 state into 
tungsten, but once emplanted, the ionization of neon was unknown. 
All combinations of W^, W ^ and W ^ with Ne^ and N ^ were sub- 
jected to Hartree-Fock analysis. (See Figure 5-) Only Ne^-W^^ 
interacted in an approximately exponential manner and could there- 
fore possess realistic BM parameters. Attempts to linearize 
o +1 o o 

Ne -W and Ne -W were made, and subsequent BM parameters were 
determined. The results of such changes did not significantly 
influence the results of this investigation. 



E. RUNNING TIME 

The following factors effected the problem running time: range 
of potential, size of crystal, depth of mobile layers, and degree 
of damping. First, the range of the potential was picked to 
include at least NN4 interactions. The range used for the tungsten 
simulation was 3-4 LU, which includes interactions out to NN4 . The 
error made by neglecting NN5 interactions was only 3 % in the 
binding energy of an interstitial, but the omission of NN5 inter- 
actions cut running time almost 105^. Second, the size of crystal 
and depth of mobile layers were picked as small as possible, for 
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reduced running time, but were at least large enough to completely 
contain the potential range. Third, the half velocity method of 
damping was used whenever possible. In general, when velocities 
were zeroed at the end of each timestep, a timestep took about ten 
seconds, and equilibrium was reached in about 300 timesteps. \Vhen 
velocities were halved, each timestep again took about ten seconds, 
but equilibrium was reached in about 150 timesteps. 

It was originally expected that increasing [dTi] would decrease 
running time. Most W-W simulation was done with [dTiI = 2%; i.e., 
the most energetic or stressed atom could travel .02 LU before the 
damping of velocities and the starting of a new timestep. When 
CdtiI was increased, the atoms moved more erratically tov/ard equi- 
librium, and vibrated there, but did not achieve equilibrium 
cignif icantly sooner. (See Figure 6.) 



F. SUMMARY 

In summary, the steps of the program are outlined: 

1. Variables are initialized, constants established, and input 
data read in. 

2. Scaling factors and time saving multipliers are calculated. 

3. Morse and BM potential functions are calculated based on 
input data. Subsequent forces, based on derivatives of these 
functions are calculated. Potential erosion and force modifi- 
cations are performed. 

4. Potential cutoff’s CroEa] , CrOEB] , CrOECJ are established 
and the smooth fitting cubic equation is placed in the g3ip. 

5. The desired crystal type, size, and orientation is built, 
and the point defect positioned (see Appendix A). 
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6. Mutual potential energies of all atoms in the crystal are 
calculated. Local potential energy is calculated. (See 
Appendix D.) 

7 . All initial positions and potential energies are printed, 
along with total potential and total kinetic energy, local 
potential energy, and the change in local potential energy. 

8. The first tiraestep is started, with an arbitrary running 

— iZi , , , 

time of 10 sec. Velocities and positions are calculated 

by the average force method, and the maximum velocity CemaxI 

and maximum force CfMAx] are found. 

9 . A new Cot], based on LENIAX] , CfMAX], and [OTl] is calculated 
for use in the next timestep. (See Appendix C.) 

10. All velocities are zeroed or halved as an energy damping 
method; and the process (8. to 10.) is repeated. 

11. At selected tiraesteps, all changes in position (CdxJ, CdyJ, 
CdzJ), velocities (Cvx], L Vy] , Cvx3), and kinetic, potential, 
and total energies [pke], CpPe] , [pTe] for each atom in the 
crystal are printed. 

12. The program is ended after a pre-selected timestep, with a 
final printout of position and potential energy of each atom, 
as in Step 7 . 
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IV. RESULTS 



KS’s experiraental data indicated that four or five interstitial 
positions in the first few layers of a tungsten lattice could be 
found that would result in different binding energies. It was soon 
found that many parameters in the program could effect the results, 
and so a systematic attempt to isolate the effects of each indivi- 
dual parameter was undertaken. 

A. THE CRYSTAL 

As explained in Appendix A, the Y = 0 plane was the crystal 
surface; the Y = 1 plane was the first layer beneath the surface, 
etc. Each atom in each layer was then designated by appropriate 
X and z coordinates. This consxruction was independent of xype 
of lattice; i.e., the surface layer in either BCC (100) or FCC 
(100) was Y = 0, etc. An interstitial that escaped the lattice 
normal to the surface travelled in a (OlO) direction, and an ion 
that escaped normal to a side travelled in either a ( lOO) or (OOl) 
direct ion . 

The dimensions of the lattice had a great bearing on the re- 
sults. In general, the larger the crystal, the more realistic the 
results, but increased computer time prevented the use of a size 
bigger than absolutely necessary. All point defects were placed 
as close to the center of each plane as possible, and the x and z 
dimensions of the lattice Cix] and [iz] were chosen to completely 
enclose a circle of radius [rOEc] from the point defect. In this 
way any point defect in the center of the lattice would not feel 
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the effect of the sides of the lattice, especially unequal numbers 
of atoms in all directions. In the y direction, the lattice was 
again built deep enough to completely contain the radius of the 
potential of a point defect placed at the center of the lattice. 

To simulate the effect of an infinitely deep lattice, the bottom 
two layers of the crystal were held immobile, but still allowed 
to interact with all mobile atoms above them. The tungsten crystal 
size used most often was a 10 x 10 x 10 cube with the bottom 
10 X 2 X 10 volume held rigid. 

Often erratic behavior in the simulation could be eliminated 
by simply increasing the crystal size. This was especially true 
for the problem of crowdion migration. 

B. CROWDION MIGRATION 

Crowdion migration is a chain reaction of single lattice ^=5te 
jumps initiated by interstitial implantation. If the chain re- 
action ends by pushing the surplus atom into an already existing 
vacancy, the interstitial-vacancy pair is called a Frenkel pair. 

A Frenkel pair can also be created dynamically by moving an atom 
from its lattice site to a nearby interstitial position, from 
where it can cause migration back to the vacancy. If the mi- 
gration cannot find a vacancy, and travels all the way to the 
surface, the surplus atom forms a "stub’\ Normally, migration 
is always in a closed packed direction; i.e., in the (ill) di- 
rection in BCC. It was discovered, however, that this rule was 
modified near a surface, since an imbalance of forces in the di- 
rection normal to the surface automatically pushed a crowdion in 
that normal direction into a stub position. In the tungsten 
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lattice, for instance, a tungsten interstitial that did not 
initiate crowdion migration would reach equilibrium in a (llO) 
split interstitial position, as previously found by EVE U?] . A 
tungsten interstitial that did initiate crowdion migration would 
sometimes migrate in a (ill) direction because of clos ed-packedness 
or sometimes in a (lOO) direction if implanted near a (100) sur- 
face. Crowdion migration was never found in a (llO) direction, 
since this is the least closed-packed of these three directions* 
(Hence the tendency toward split-interstitials in this 
dir ect ion ) , 

Crowdion migration was a very common process near a lattice 
surface. It was found, however, that varying the choice of atoms, 
the range of the potential, and the rate of energy damping could 



enhance or reduce the tendency tcivard crcv;dicn mi< 
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already mentioned was the fact that increased crystal size re- 
reduced crowdion migration. Particular attention was paid to the 
proper choice of values for these parameters, in order to cor- 
rectly determine whether or not crowdion migration actually existed 
This question of crowdion migration was especially critical in this 
simulation, since the binding energy of a particular atom is a di- 
rect function of the nearness of * its neighbors. An atom in a 
split interstitial position feels a more repulsive potential than 
an interstitial that has initiated crowdion migration and 
’stolen’ a lattice site and has thus reformed the original perfect 
lattice with every atom in a normal lattice site. 

1 . Choice of Atom 

Some elements tended to initiate crowdion migration more 
than others. This applied to both choice of lattice atom, and 
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choice of interstitial atom. For instance, crowdion migration 
was much more common in tungsten than in copper. This was due 
to both the size of the tungsten atom and the nature of the crystal. 
Also it was found that different element point defects in the same 
lattice produced varying degrees of crowdion migration. A neon 
interstitial is so small and light that it never initiated crowd- 
ion migration, even when only one layer separated it from the sur- 
face. An argon interstitial initiated crowdion migration at the 
surface, but not deeper in the lattice. A tungsten interstitial 
always initiated crowdion migration, unless placed at the center 
of a huge l4 x l4 x l4 tungsten lattice. This is an example of 
increasing crystal size to prevent crowdion migration. In 
general it can be stated: the more massive the interstitial, the 
mbre probable crowdion migration. 

2 . Range of Potential 

In surface simulation, the range of the potential was a 
critical factor, since it determined whether or not an atom could 
"see” the surface. Because copper has been the standard element 
for lattice simulations, many versions of a copper potential with 
various ranges, have been determined. As mentioned previously, 

GW C 12 ] calculated a Morse potential for copper that effectively 
had an infinite range (150th nearest neighbor). If this potential 
is truncated at very close ranges, i.e,, NNl or NN2 , the potential 
is seriously underestimated. This under estimate rapidly dimin- 
ishes as the truncation range increases. Since GW parameters for 
the Morse potential could not be used for NN2 interactions, 

Anderman Cl4l calculated parameters for a Morse copper potential 
that would approximate GW^s results in simulations truncated after 
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NN 2 • He did this by deepening and broadening the well. Although 
Anderman parameters and GW parameters led to very similar results 
in an infinite lattice, they led to quite different results in this 
simulation. In general, if the range of a point defect potential 
function overlapped a surface, crowdion migration would take place 
toward that surface, because of an imbalance of forces in the nor- 
mal direction. The effect was a little more complex than this be- 
cause of surface relaxation: if the range of the point defect po- 
tential function overlapped a relaxed surface layer, the slight 
force imbalance would again result in crowdion migration. Accord- 
ing to GGMV, "the machine calculation showed that this atom rapidly 
moved... in a direction determined by minor asymmetries in the 
starting conditions...**^. In this copper simulation, an inter- 
stitial placed in the forth layer with a GW potential range of 
3.1 LU (NN 4 ) caused complete crowdion migration, resulting in a 
copper stub on the surface. An identical run with Anderman para- 
meters for an NN 2 potential to a range of 2.4 LU resulted in a 
(100) split interstitial with minor, damped migration to the sur- 
face. Instead of a stub copper atom as before, four copper atoms 
in the surface layer bulged about .4 LU. 

Another example of a short range potential which demon- 
strated this lack of ability to initiate crowdion migration was 
the repulsive foreign defect potential, with a range of \P 3 " LU. 

In copper, this potential quickly led to split interstitial 
positions and no crowdion migration for all copper interstitials 

^Gibson, J.B., Goland, A.N., Milgram, M., and Vineyard, E.H., 
'^Dynamics of Radiation Damage, * *The Physical Review , V. 120 , Xo. 4 ; 
p 1237, Nov 15, i960. 
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except those placed in the first two layers. In tungsten, even 
this short range potential could not retard crowdion migration in 
the 10 X 10 X 10 lattice. Only in the center of a l4 x l4 x l4 was 
a tungsten split interstitial stable. This stability applied only 
to the short range repulsive potential: a repeat run using the 
standard composite potential with a range of 3.4 LU initiated 
crowdion migration. This increased range enabled the interstitial 
to find minor asymmetries in even a l4 x l4 x l4 lattice. 

3 . Energy Damping 

Energy damping was accomplished in this simulation by 
reducing each atom's velocity at the end of each timestep. Two 
methods were used: at first, each velocity component of every 
atom in the crystal was zeroed at the end of each timestep. Later, 
the halving of each velocity component at the end of each time- 
step was employed to save computer time. In a tungsten lattice 
the results of both methods were the same; all final positions and 
binding energies were identical. Neither method prevented crowd- 
ion migration. In copper, these two methods led to si ight ly ' dif- 
ferent results. Although the final position and binding energy 
of an interstitial was almost identical, and although crowdion 
migration was initiated in both cases (GW's parameters and a 
3*1 LU range were used), the zeroed velocity method had damped the 
migration significantly by the time it reached the surface, where- 
as the halved velocity method caused a complete, undamped mi- 
gration to the surface. 
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C. INTERSTITIAL IMPLANTATION 



All interstitials were placed in the obvious holes in a hard 
sphere, close-packed lattice model. (See Figure 7 .) Every 
"hole” in the tungsten lattice had exactly the same geometry; 
i.e., two neighbors 1 LU away; four neighbors 2 LU away, four 
neighbors 3 LU away, etc. The only factors which differentiated 
between these identical holes and thus led to different binding 
energies were: layer number, or lattice depth, and open channel 
direction. An interstitial in the third layer was more tightly 
bound than one in layer two, etc. Also, an atom in a given layer 
could be placed in two types of holes: one in which the inter- 
stitial was in the BCC (OlO) open channel direction, in which case 
the interstitial could ”see” the surface; and one in which the 
interstitial was in the BCC ( lOO) or (OOl) open channel, in which 
the interstitial could not "see” the surface. Note, however, that 
if an atom could not ”see” the surface, then there was no difference 
between these two positions, since both have two neighbors 1 LU 
away, four neighbors 2 LU away, etc. Note also that even if these 
two sites are identical and possess exactly the same binding 
energies, the difference might still show up in diffusion pro- 
bababilities : an interstitial in a (OlO) open channel in the 
second layer must only move two lattice units to escape the cry- 
stal. An interstitial in a <100) or (OOl) open channel must move 

1 LU in either the x or z direction into an open channel, and then 

2 LU to escape; i.e., it must move like a knight in chess. This 
extra step might lead to a different diffusion probability. 
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D. THE TUNGSTEN LATTICE SELF DEFECT 



The tungsten self interstitials and self replacement defects 
were the chosen standard for this analysis, as explained in the 
in the OBJECTIVE . The tungsten defects could be treated as 
lattice atoms and allowed to interact with all other atoms with 
the composite potential (Method 1); or they could be treated as 
foreign defects and allowed to interact with only a repulsive 
potential (Method 2). A total of three different defect positions 
were simulated: an interstitial in a (OlO) open channel ( int A), 
an interstitial in a (lOO) or (oOl) open channel (int B), and a 
replacement atom (rep) in a lattice site. (See Figure 7-) 

1 . Interstitials 

As previously mentioned, all tungsten interstitials initi- 
ated crowdion migration v;hcn treated by method 1. An interstitial 
treated by method 2 also initiated crowdion migration unless buried 
in the center of an enlarged l4 x l4 x l4 lattice. Because an 
interstitial that has pushed its neighbors away has a lower po- 
tential than one that has not done so, the numerical values for 
binding energy of tungsten interstitials could not serve as a 
true standard for comparison with W-Ne results. Qualitatively, 
however, much could be learned from the W-W energy levels. First, 
it was expected that all energy levels of a defect found by 
Method 1 would be negative at equilibrium. Values for the compo- 
site potential can be either negative or positive, but are posi- 
tive only at very small separations. A negative potential energy 
means the atom is bound in the crystal. As shown in Figure 1, on 
the next page, the first two interstitial levels for the V7-W 
reaction. Method 1, were at -3-0 eV and -3.1 eV, corresponding to 
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the two interstitials in the first layer (int A and int B) . The 
next two levels were at -5.4 eV and -5.9 eV, corresponding to the 
interstitials in the second layer. Interstitials in the third 
layer and deeper had binding energies ranging from -7.4 eV to 
-8.1 eV. The value of -8.1 eV, labeled was obtained from the 

interstitial placed in the center of the seventh layer of a 
l4 X l4 X l4 lattice. Since this value, and all other values for 
Method 1 binding energies were reached after crowdion migration, 
they were all expected to be lower than they would be without 
crowdion migration. The only way the true binding energy of a 
tungsten interstitial could have been found would have been to 
find a tungsten crystal size large enough to contain the crowdion 
migration of a tungsten interstitial with the long-range composite 
potential. Compute running time made this impossible. 

Also shown in Figure 1 are the binding energies for the 
Method 2 W-W interstitials. Note, first, that they were all 
positive. This was again expected, since the potential equation 
for Method 2 is positive over all space. Note, second, that the 
binding energies for interstitials in the first two layers were 
zero. This was because all purely repulsive atoms in these first 
layers escaped the lattice completely. The other positive levels 
shown, were incomplete, because many interstitial positions did not 
possess stable energy levels. The unstable levels oscillated be- 
cause of significant lattice motion, caused by crowdion migration, 
and measured by a short range potential. The range was so short, 
that significant jumps in binding energy occured when an atom 
moved into, or out of the range of the potential. Evidently, 
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small vibrations in atoms with an equilibrium distance of about 
^3" LU from the interstitial, frequently caused crossings of this 
range limit, adding or subtracting energy from the binding energy 
each time one of them crossed, thus invalidating many of the inter- 
stitial binding energies. 

The energy levels shown, however, demonstrated the meaning 
of a positive ^^binding energy". The numbers reflected the "amount 
of repulsion" associated with different positions in the lattice. 
The ordering of the levels, i.e., higher energies for deeper layers 
was expected. An interstitial deeper in the lattice felt more 
repulsion, because it was surrounded by a greater number of re- 
pulsive neighbors. Again, the level labeled "«" represented an 
interstitial in a l4 x l4 x l4 lattice; but in this case of a 
Method 2 inter s+ i+ ia 1 . crowdion migration did not occur. 

The concept of a positive binding energy may or may not 
be the actual physical situation, but it is still academically 
valuable. Instead of an atom resting near the bottom of a po- 
tential well, as in Method 1, an atom can be "wedged" between the 
repulsive walls of its neighbors. In both cases, the atom is 
"bound" in the lattice. The ordering, spacing, and other cor- 
respondences between the positive and negative levels validate the 
qualitative use of the positive levels. 

2 . Replacement Impur it ies 

A tungsten lattice with a replacement impurity is a per- 
fect tungsten lattice, but Method 1 or Method 2 could be used on 
the atom. in question. For Method 1, a perfect crystal wasallowed 
to relax with time, yielding only negligible motion; and the 
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binding energies of the center atom in each layer was recorded 
(see Figure 1). The binding energies of atoms in the first, 
second, and third layers were-5.3 eV,“ 6.8 eV, and-7.0 eV respec- 
tively. The subsequent reversal of order for levels corresponding 
to deeper layers is a program anomally, caused by the use of a 
finite depth crystal. Runs on larger crystals indicated that the 
order v;ould not reverse in an infinite lattice. The level 

at-^ 8.8 eV is the experimentally determined heat of sublimation [ 26 ] 
of tungsten. These numerical values for the binding energies 

were valid as standards of comparison for the Ne-W data, since the 

motion and equilibrium positions of the replacement atoms and their 
neighbors were nearly identical, and usually less than .3 LU in 
both cases. Note that the binding energies of the Method 1 re- 
placement atoms were lower than the interstitial atom levels. As 
previously stated, this was to be expected, since a replacement 
atom rests in the bottom of a periodic potential well in the 
lattice, whereas an interstitial rests in a higher well, because 
it is nearer to its neighbors than the normal equilibrium sepa- 
ration. Also note that if the entire Method 1 spectrum of binding 

energy levels were used as a standard of comparison for Ne-W levels, 

then the W-W interstitial levels .should be higher with respect to 
the replacement levels, than shown in Figure 1, because of crowdion 
migration . 

The Method 2 replacement level labeled cor responded 

to a replacement atom in the center of the fifth layer in a 10 x 
10 X 10 lattice. Note that it was not above the interstitial 
levels, as it should have been by comparison with Method 1. This 



35 



was the major weakness of Method 2: because it was a measure of 
repulsion only, and neglected the potential well, it under esti- 
mated the values of the binding energies of atoms whose normal 
position was in that well; i.e., replacement defects. 

E. THE NEON DEFECT IN TUNGSTEN 

The neon atom defect was again placed in any one of the three 
lattice positions, labeled ”int ”int B** , and ”rep”. The neon 

defect could be treated by Method 2 only. The neon energy levels 
are shown in Figure 2, on the next page. Note that the energies 
have been multiplied by a mass correction factor. (See Appendix D. ) 

1 . Inter st itials 

The neon interstitial never initiated crowdion migration, 
but as in the W-W case, atoms placed in the first two layers es- 
caped the crystal, and therefore had zero binding energy. Again, 
the ordering of the levels was a measure of the replusion on the 
interstitial, which increased as the interstitial was placed 
deeper in the lattice. Again the level was the result of 

placing a neon interstitial in the center of a l4 x l4 x l4 tung- 
sten lattice. 

2 . Replacement impurities 

Again note that the replacement levels are lower than 
would be expected by comparison to the W-W Method 1 standard. It 
is hypothisized that these replacement levels should be higher 
than the interstitial levels, by comparison to the standard. This 
assumption is valid if a Ne-W potential well exists. It is fea- 
sible that since neon is almost incapable of binding, that no Ne-V.' 
well exists. On the other hand, the ionization state of both neon 
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TEMPERATURE, IOC 



and tungsten in the lattice is unknown, so the existence of a 
shallow potential well is quite possible. 

F. CORRELATION WITH EXPERIMENT 

1 . Scaling the Levels and Peaks 

since no direct way of transforming the Ne-W energy levels 
into correctly scaled negative values exists, an arbitrary linear 
scaling factor between KS^s data and the Ne-W levels has been 
used. Note that every level or group of levels corresponded to 
an experimental peak in Figure 2, except in two places: first, 

the broad peak at about 2000^ K had no energy level counterpart, 
but could be assum.ed to correspond to the closely ordered replace- 
ment levels, shifted above the interstitial levels by the above 
hypothesis. Second, the narrow peak at 450 was without an 
energy level counterpart. Note that three interstitial levels 
had zero simulated binding energy because they had escaped the 
crystal. If, however, the assumption that a Ne-W well exists was 
true, then some or all of these three interstitial locations, 
(i.e., two surface layer positions and the open channel position 
in the second layer) would be stable, bound positions, and would 
be expected to generate an energy level in the vacinity of the 
450^ K peak. If the arbitrary scaling of the peaks and levels was 
correct, then the peak at 450^ K is proof that a Ne-W well exists, 
since a purely repulsive potential would not allow a first layer 
or second layer open channel interstitial energy level to exist. 
The existence of this well, then, would in turn substantiate the 
shift of the replaceraent levels above the interstitial levels, to 
correspond to the 2000^ K peak. 
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Probes of Potential Wells 



Various attempts to substantiate the scaling between the 
levels and peaks were made. No approximation to a Ne-W potential 
well could be justified, and the correspondence between W-W and 
Ne-W results was not complete enough to invert and scale the po* 
tential energy levels to realistic, negative, binding energies. 

Attempts were also made to match both Method 1 and Method 2 
W-W energy levels to KS’s data. The Method 2 levels were too 
incomplete, and the Method 1 levels required an unknown arbitrary 
reduction of interstitial energy levels to compensate for crowdion 
migration. Too many alternate reductions were possible to choose 
one as the correct rescaling. 

Another approach that yielded little information was a 
plot of the differences belweeii iiiterst itial arjd replacement 
energy levels for each layer in the lattice. 

One valuable method of investigating the nature of possible 
Ne-W negative binding energies was to probe the perfect lattice 
with an interstitial at various initial positions and plot the re- 
sultant potential energy of an interstitial vs. position. Since 
the potential was always positive, the results of this investigation 
were potential wells above the x-axis. Although the positive lo- 
cation of these wells was not realistic, the relative depth of the 
wells was significant. The average depth of a neon interstitial 
well, deep in the lattice, was about 4.2 eV (see Figure 8). The 
graph was made by placing interstitials in positions in the (OlO) 
open channel, and thus represents the barriers that the interstitial 
must penetrate as it escapes the crystal. Note that the w 11s were 
not at the obvious holes in the BCC lattice, but were between layers. 
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In the actual simulation, the interstitial rarely fell into this 
well, but instead pushed its two nearest neighbors away and made 

initial position the low potential position. Here again, sur- 
face effects and relaxation reduced the tendency for interstitials 
to relax into expected infinite lattice equilibrium positions. 
Slight differences in the final equilibrium positions were not of 
significant importance to binding energies. The actual numerical 
values for the depths of these positive wells were not necessarily 
scaled properly since they ignored the actual Ne-W potential well, 
and were found from perfect lattice probes at time zero, before 
relaxation. Nevertheless, the depth of a well deep in the lattice 
of 4.2 eV agreed well with KS’s prediction of 4.5 eV [8] for the 
desorbtion energy corresponding to 1720^K, at the front edge of 
the highest peak. 1720^ K closely approximated the temperature 
that the arbitrary scaling had assigned to a deep interstitial. 

Note that the initial position of a replacement atom was 
\[^ LU from its neighbors, and thus because of potential erosion, 
it had a potential of zero eV initially. To climb out of this 
well, about 17.5 eV must be supplied. Although this number was 
inaccurate for the same reasons listed above, it did demonstrate 
that even a purely repulsive potential can predict a greater 
binding energy for replacement atoms than for interstitial atoms. 



V. CONCLUSIONS 



An arbitrary scaling has been used to correlate the simulation 
results with experimental data. Although the method was not ana- 
lytically sound, no other avenues of approach to the problem could 
be found that could further justify our hypothesis. 

Satisfaction can be gained, however, from the fact that these 
results compare favorably with known data at many interfaces. Our 
model was a tried and proven one, with many successful sputtering, 
channelling, and similar simulations to its credit. This present 
model invariably behaved in a physically valid manner or a manner 
which could be made physically acceptable by varying the con- 
trolling parameters in the program. Specifically, many previous 
experimental and simulated results for infinite crystals were 
reproduced when simulation took place deep in a large lattice, 
such as the (lio) split interstitial position for BCC structures. 
The Method 1 replacement levels, if found for a much deeper cry- 
stal, would have asymptotically approached very close to the 
8,8 eV heat of sublimation. The simulated depth of the positive 
potential well for interstitials of about 4*2 eV closely approxi- 
mated KS*s prediction of 4.5 eV . 

All avenues in additional computer simulation have not been 
exhausted. Future simulation of Argon, Krypton, and Xenon defects 
in tungsten should be fruitful. Comparisons between the relative 
locations of these new energy levels might further substantiate 
this research. In particular, if simulation can explain why KS’s 
neon data contains five desorbtion peaks, while their Argon, 
Krypton, and Xenon data contain four peaks, it will be a major 



success. KS also gathered data for different crystal surfaces 
and different angles of incidence, which might be investigated 
by computer simulation. 

This simulation was also important in that it investigated 
the lattice surface; a topic which has not received as much 
attention as infinite crystal dynamics. Since radiation damage 
theory and modern transistor theory is very much concerned with 
the crystal surface, the computer simulation field vjill undoubtably 
increase their emphasis on surface effects, with considerable 
attention toward better ways of treating a foreign interstitial. 

A new exhaustive book which reports on the present state of 
knowledge in all these areas, with emphasis on experimental re- 
sults has gust been published. It is a report on the proceedings 
of the International Conference on Vacancies and Interstitials in 
Metals, 1968 [27] . 
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APPENDIX A: CRYSTAL GEOMETRY 



The computer program can call any one of nine lattice generator 
subroutines: three face-centered cubic subroutines ([lIOO] , ClIIO], 
ClIII]), three body-centered cubic subroutines ([bIOO], CbIIO], 
[Bill]), and three diamond subroutines ([dIOO], CdIIO], [d111]). 

The diamond subroutines are never used and therefore not compiled; 
but provision has been made for their future inclusion in the pro- 
gram. The dimensions of the lattice chosen were controlled by the 
input data variables [ix], Liy], and C iz] . Each atom in the cry- 
stal was numbered, in the order x followed by 2 , followed by y. 

For the surface layer (Y 0) of the tungsten 10 x 10 x 10 lattice, 
atoms were numbered from 2 - 26 ; for the first layer below the sur- 
face, atoms were numbered 27-51? etc. Atom number 1 was the pri- 
mary, 01 point defect atom. Lld] was the number of the last mobile 
atom, or the last atom in the eighth layer, number 200; and [ll 3 
was the number of the last atom in the crystal, number 250 . 

The placement of point defects was accomplished as follows: 

after the desired perfect lattice was built to the desired size, 
subroutine PLACE was called. Three types of defects were allowed: 
vacancies, interstitials, and replacement impurities. The type 
and location of the defect were controlled by input data vari- 
ables: the type of defect [iTYPe], an atom number [nVAC] and a 

displacement vector CdIX, DlY, Dlzl in LU . If CiTYPEJ = 1, a 

vacancy was created in site number CnVAC] . This "removal" was 
accomplished by setting ClCUT (NVAC)] = 1 which "turned off" the 
atom, removing it from all calculations. If CiTYPE] == 2, an 
interstitial was created in a position C-DIX, -DIY, + DIZj LU from 
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site number f.NVAC]. This interstitial was always atom number 1. 
Since atom number 2 was alivays at the origin, number 1 could be 
placed using a displacement Vector from the origin (from [nvac] = 2) 
or using a displacement vector from a site next to the interstitial. 
If L ITYPe] = 3j 2k vacancy was created in site number CnVAC] and 
a replacement impurity, put in its place. Note that for both 
[iTYPEJ = 2 and 3? either a foreign or self defect could be 
placed. For the case of either the self-interstitial or the self- 
replacement atom (giving us back the perfect crystal), either 
method 1 or method 2 of calculating the potential could be used. 

The choice of methods was also an input parameter: for method 1, 

Ciq] = 1, and for method 2, [ IQ] = 2. 



44 



APPENDIX B: POTENTIALS 



(In this appendix and all subsequent appendices, the brackets 
denoting program language are dropped. Program language is still 
written in all capital letters.) 

A. BORN-MAYER REPULSIVE POTENTIAL 

1. Potential Energy : For the lattice atom interactions, the 
Born Mayer potential equation is; 



goes to zero at the nearest neighbor distance. in the program. 




( 1 ) 



or 



POT = EXP(EXA. + EXE-X'DIST) 



(lA) 



For bullet-lattice atom interactions, 




(3) 





POT = EXP(PEXA + PEXB*DIST) - PPTC 



(3A) 



2. Force: For the lattice atom interaction, 




= -B exp (A + Br ) 




in the program, Ct2(-B + A) = ( ALOgC -EXB*CVEd] + EXA) = FXA, 
where CVED is a conversion factor for units, and 



FORCE = CfXA + EXB*DIST] 



(^A) 



For bullet-lattice atom forces. 



= 5r . . = 





where (2n-B ' +A ' ) =( ALOgC -PEXB*CVED] +PEXA ) =PFXA, and 
FORCE = EXpCpFXA + PEXB*DIStJ 



(5A) 
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f 

c 






for buUe.-lattrce ato™ forces for K- 

ces for which lattiro 
range of the bullet's f ' 

force during the tinestep a fi 

mation to the forc^a -ic • approxi- 

ice IS given by.- 



Force = 

^ij~ROE 



(6) 



(See Figure o ) v w 

ah * retarded potential • 

Oboue. . conuorsion fagtor is „oeded to preserue th 

- ^ohstitutions are pade in the ret 

— " <pe., hecooes <. „eOtp.», , pi 

Vf, = EXP[pactPHXB.Msx3, and V .(POE, = Ppfc p ’ 

- Vt^.iPOE, calculated after th Pac p 

'if- = 013X.P0E = 1.P, and 



force 



= - ppprr 

DFF ^ 



(6a) 



N°te: ( 5 A) is used for 0 ^ prs^ ^ ROE - OTl - r 
( A) IS used for RoE^, < ^ 

’ * r • ^ 

MORSE POTENTIAL 

-£^i2ntiaI__Enerqy; p^r 

potential equations is, ■ '"*^‘'‘''=''ons , the Morse 

ij “‘■®'^f'^(‘^ij-r^)).2epp{^(r.^.-rJJ (2, 

= exp[ (*!D + 2ttr i./sq.,, 7 , 

o' <-»)tijJ-expC(fc( 2 D)t Ctt^).(a,r 3 ■■, - - ‘ ' 

= EXpf(,d,OG(OCON)t ° . • ■ V .: V ;• 'V 

- EXpC(AL0Gt2 .nc ALPHA*CVR) .DISt] 

= expCcgdi C ™’"' '^™''*EE) -(alpha..cvr).distJ- 

LCGD1-CGB1*Dtct1 r ■ 

distJ- expLcgd2-cgb2*dist] 

2 • (-2A) 

• atoM interactions/ , 

Force = ^ 1 «. nT‘ 

6r ^20/ exp[-2a(r . ,-r ) ‘I 2r, f 

ij o^-^ 2a expi-a(r. -r )'■ 

ij o' 






.* 



^6 



I 





= EXpCaLOG( 2*ALPHA*CVR*CVED) +AL0G( DCON) +2 . *ALPHA*RE 

- ( 2 , *ALPHA*CVR ) *DISx] 

- EXP[ ALOG( ALPHA*CVR*CVED) +ALOG( 2 . * BOON )+ ALP HA* RE 

- ( ALPHA*CVR ) *DIST'J 

= EXP[ ALOG( -CGB1*CVED) +CGD1 ) -CGBl*DISx] 

- EXpL' (AL0G( -CGB2*CVED)+CGD2)-CGB2*DISX] 

= EXP[CGF1-CGB1*DISXJ- EXP[CGF2 -CGB2*DISx] (7A) 

C. CUBIC FIX 

1. Potential Energy : The best cubic fit between the BM and 
Morse potentials is calculated in Subroutine CROSYM* The po- 
tential equation, defined between ROEA and ROEB, is 



POT = CP3r . r + CP2r . . 



+ CPlr . . - 1 - CP0 
j-J 



or 



POX + DISX*(DISX*(DISX*CP3+CP2)+CP1)+CP0 (8a) 



2. Force: Force 



-5POX 
br . . 




CPI (9) 




ij ij 



FORCE = DISX*(DISX*CF2+CF1)+CF0. 



(9A) 
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APPENDIX C: AVERAGE FORCE METHOD AND TIME DURATION THEORY 



A- AVERAGE FORCE METHOD: The average force technique has been 
explained in great detail in Ref. 15. It was summarized in 
Chapter 3> smd therefore discussion here is limited to the average 
force method in the program language. 

When the desired lattice is built, the position of the i1U|i atom is 
stored simultaneously in RX(I), RY(I), RZ(I); RXK(I), RYK(I), 
R2K(I); and RXI(I), RYI(I), RZI( I) . The latter set of coordinates 
never change and are used for comparing new positions to original 
positions, and in calculating DX(I), DY(I) and DZ(I) for output. 
The middle set of coordinates containing the letter K are for 
storing the initial positions at the beginning of each timestep. 

A step by step summary of the average force method, showing X 
coordinate calculations only, follows: 

1. Based on the position, RX(I), of the i th particle at the be- 
ginning of the timestep, the force FX(I) is calculated in STEP. 

2. RX(I) is stored in RXIC(I). 

3. The new, temporary position, RX(I) is calculated, based on 
the force at RXK(I): 

= X + VAt + F(At)^/2m (10) 

or 

RX(I) = RX( I)+DT0D>'(HDT0M*FX( T)+VX( I) . (lOA) 

4. A new force is calculated, based on this new position RX(I). 

5. VSS stores VX(I), the or iqinal velocity of the i th atom. This 
velocity is half the velocity of the i th atom in the previous 
timestep : the H factor being an arbitrary damping multiplier. A 
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( 11 ) 



new velocity, based on the new force, is found: 

= V + FAt/2m 



or 

VX(I) = VSS+HDTOM*FX( I) . (HA) 

6. The final position is calculated, based on the average of 
these two velocities 

X^ = X + %At(V+V^) (12) 

RX(I) = RXK( I) +(VX( I)+VSS)*HDTOD (12A) 

The resultant velocities are halved, a new timestdp duration is 
calculated, and the process repeated. 

B. TIMESTEP IXJRATION THEORY 



This simulation uses the best possible estimate of a timestep 
duration, DT, as calculated from the present s tate of the forces 
and energies in the lattice for use in the next timestep. To 
limit motion to an increment small enough to preserve the accuracy 
of the average force approximation, we define DTI as the maximum 
distance any atom is allowed to move in one timestep. 

From (10), we find 



Therefore, 



For 



AX. = (V.+ F.At/2m)At. 
1^11 ' 

At = AX./(V.+F. At/2m) . 

i'll ' 

V. » F.At/2m, At = AX./V.. 
1 1 ’ 11 



(13) 

( 14 ) 



If we find the fastest moving atom and assure that it does not 
move more than DII, we have limited the motion of all other atoms 
to less than DTI. 

DT = (DTI*CVD)/EMAX = FDTT/EMAX (14a) 
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m 





where 

EMAX = SQRT(VX(I)*VS(I)+VY(I)*VY(I)+VZ(I)*VZ(I) . 

For » F^At/2m, 

2mAX. 

At = AX^/F^At/2m = — . (15) 

i 

Anatagous to above, we find the most stressed atom and assure 
^that it does not move more than DTI. Thus, 

DT = SQRtC (2.*PTMAS*DTI*CVI))/FMAX] 

= sqrt[tfac/fmax] (15A) 

Where FMAX= SQRt[ FX( I ) *FX( I ) +FY( I) *FY( I) +FZ( I ) *FZ( I ) ‘J . 

Since rigorously we cannot make either or these limiting assumpt- 
ions, we must go back to our original equation for DT, equation (13)* 
Since this equation involves DT , we proceed as follows; 

1. Assume « F^At/2m and calculate At from ( 15 a) . 

2. Insert this preliminary value for DT in (13) compare 

Vf to F^At/2m. If is larger, calculcxte DT from (l^A). If 
F^At/2m is larger calculate DT from (15A). 

A complication arises when a foreign impurity is in the lattice, 
because of a variation in the value for m in (15). This is 
especially acute when the differences in masses are great. The 
method used to solve this problem is as follows: 

1. If eithe r FMAx = ££ EMA X = , the entire proceedure, above, 

is followed using the mass of the bullet for m. 

2. If both FMAX 7^ F^ and EMAX f , the entire proceedure is fol- 
lowed using the mass of a lattice atom. 

The requirement that the bullet mass be used if either EMAX or 
FMAX describe the bullet circumvents the problem of having the 
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bullet the fastest moving atom, but not the most stressed atom, 
or visa versa. 

. . ” l4 

To begin the problem, an arbitrary value of 10 seconds is 
assigned to DT. If at any time in the program EMAX = FMAX = zero, 

-l4 ... ... 

10 IS again assigned to DT to prevent division by zero. 
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APPENDIX D: SUBROUTINES STEP, ENERGY, AND LOCPX 



A. DISTANCE CALCULATIONS 

In all three subroutines, STEP, ENERGY, and LOCAL, a method of 
finding all atoms within a given radius of another atom was needed. 
For lattice atom interactions, atoms inside ROEA, ROEB, and ROEC 
were found; for the foreign interstitial interactions, atoms inside 
ROE were found; and for LOCAL, atoms inside ROEL of a point defect 
were found. The t ime saving technique used to do this was to 
successively eliminate all atoms with an x component difference 
greater than the given radius, then similarly for y components, 
then for 2 components. The resulting volume not eliminated is a 
cube circumscribing the desired sphere. Finally the time-con- 
suming test of eliminating all atoms for which the desired radius 
is less than SQRT (DRX^DRX DRV^DRY s- DRZ^DRZ) is applied to 
only atoms inside the cube. 

B. SUMMATION INDICES IP AND IQ 

Interactions V. 0. and F. . are found by evaluating all values 
ij 11 11 



in the half matrix. For example ’ ^13’ ^l4 



are found, then 



F 23 J ^ 34 ’ etc. The variable IP controls the starting point 

for the j summation. IP is always set to I + 1 to avoid the re- 
petition of finding F^^ and iQ controls the starting point 

for the i summation. If the primary is to be treated as a lattice 
atom, IQ = 1. If it is to be treated as a foreign particle, 

IQ = 2, and all F^^ are found separately. 



C. DISTRIBUTION OF FORCES AND POTENTIAL ENERGIES 

The forces F^^ are equal and opposite on i and j: i.e., * 

The potential energies are split in proportion to the reduced 
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mass of the interacting particles. This is easily understood by 
observing the kinetic energy distribution of an elastic collision 
of m and M where M ^ m. We find that m carries away almost all 

/ M \ 

the kinetic energy: specifically it carries away i — 77 

\ m+M / lOr lOL 

If a pair of atoms are to behave elastically, the potential 
energies which are transformed into kinetic energies of motion 
must be split in the same manner. For this reason, 

( Vm lr f r A ~ ) * = BSAVE.POI; 

and 

= TSAVE.POT. 

note, for BMAS = TMAS, BSAVE = TSAVE = and the energies are 
split equally. 

D. GUBI^OUTINE LOCAL 

LOCAL measures the change in potential energy associated with a 
sphere of radius ROEL surrounding a point defect. It sums up the 
potential energies of each atom found inside this^here at time 
zero. It remembers these atoms, and for each timestep re-suras 
the potential energies of these same atoms. The sum, total local 
potential energy TLPE, is subtracted from TLPE at time zero 
(TLPE 0 ), to give a measure of the change in potential energy 
(DLPE) inside ROEL. 
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APPENDIX E: COMPUTER PROGRAM GLOSSARY 



NOTE: In this glossary, the terms "point defect atom", "bullet", 

and "primary" are synonymous; and the terms "lattice atom" and 
"target" are synonymous. 

ALPHA: Input Morse potential parameter 

BSAVE: Target mass/( target mass + bullet mass); distributes 

potential energy between target and bullet 
BIND: Negative of the total potential energy (TPOT) at time zero 

BMAS: Mass of bullet in amu 

BULLET: AlphcX -numer ic array for point defect material 
CFO, CFl j CF2 : • Force parameters of cubic fit between Morse and 
Bor n-May er functions 
CGBl 5 CGB2: Morse potential parameters 

CGDl 5 CGD2: Morse potential parameters 

CGFl , CGF2: Morse force parameters 

CPO, CPI, CP2, CP3: Potential parameters of cubic fit between 
Morse and Born-Mayer functions 



CVD: 


CVR 


X 10 “^°, 


converts 


lattice units 


to meters 


CVE: 


1.6 


X 


con ver ts 


electron volts 


to joules 



CVED: CVE/CVD, a ratio used to avoid repeated division 
“27 

CVM: 1.672 X 10 , converts atomic mass units to kilograms 

CVR : LU in angstroms; converts lattice units to angstrom units 
DlX, DlY, DlZ: Displacement coordinates for location of interstitial 
from reference atom , NVAC 
DCON: Input Morse potential parameter 

DFF: ROE-DIST, the distance closer than ROE that an atom is to 

the primary 
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DIST: Distance between any two atoms 

IX.PE: TLPE-TLPE0, the change in total local potential energy 

since time zero 

DRX, DRY, DR2: x,y,z components of DIST 
DT: Length of a timestep in seconds 

DTI: Number of lattice units most energetic atom may move in 

one timestep 

DTOD: DT/CVD~-a ratio used to avoid repeated division 

DTOM: DT/PTMAS--a ratio used to avoid repeated division 

DTOMB: DT/PEMAS--a ratio used to avoid repeated division 

DX(I), DY(I), DZ(I): Change in position of iUi atom from initial 
position at time zero 

.EMAX: The maximum energy encountered in any cycle 

EV; Primary energy in electron volts 

EVR : Primary energy in kilo-electron volts 

EXA, EXB: Input Born-Mayer potential function parameters for the 
tar get 

F2 : Square of theibrce on a specific atom 

FA: The component force increment on an atom 

FDTI: DTI X CVD, a parameter used to determine DT ray maximum 
energy method 

FN: A small number used in checking potential energy zero 

point 

FM2: FM squared 

FMAX: Maximum total force on the most stressed atom in the 

cr ys tal 

FOD; FORCE/DIST--a ratio used to avoid repeated division 
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FORCE: Numerical value of the force function with a variable 



parameter 



FX(I), 


FY(I), FZ(I): x,y 


,z components 


oi' total 


force on an 


at ora 


FXA: 


Born-Mayer force 


function parameter 






HBMAS : 


% BMAS--a 


ratio 


used to avoid 


rep ea ted 


divis ion 




HDTOD; 


% DTOD--a 


ratio 


used to avoid 


repeated 


divis ion 




HDTOM : 


% DTOM--a 


ratio 


used to avoid 


repeated 


divis ion 




HDTOMB: 


^ DTOMB-- 


a ratio 


used to avoid 


repeated division 




HTMAS : 


% TMAS--a 


ratio 


used to avoid 


repeated 


divis ion 




11: 


Var iable 


in cubic fit subroutine 






13: 


1! 


ti It 


It tt 








I DEEP: 


Number of 


mobile 


layers 








IHl: 


Alpha numeric array for program title 






IH2 : 


!! t! 


1 1 


” Morse 


f unct icr* 


w ^ 




IHB: 


11 n 


tt 


” bullet 


element 






IHS: 


ft t! 


It 


” type and orientation of crystal 


IHT: 


fl ft 


ft 


" target 


element 






I LAY: 


Same as IDEEP 










IN: 


Odd-even 


integer 


used to determine atom site establishment 


IP: 


Subscr ipt 


value 


of atom. Used in subroutines STEP 


and 



ENERGY 

IQ: Parameter that determines whether or not a self defect is 

to be given a repulsive potential or a composite attractive- 
repulsive potential 

ISHUT: A parameter used to shut down the program 

IT: Unsealed fixed point x coordinate used in lattice generation 

ITT: Odd-even integer used to determine atom site estab^ ishr.ent 
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ITYPE: 



ITYPE: 


Parameter used to determine the type of point defect: 
vacancy, interstitial, or replacement 


IX, lY, 


12: Number of x,y,z planes of crystal 


to 


Variable in the cubic fit subroutine 



JJ: Parameter in the BCC(lll) lattice generation subroutine 



JT: 


Unsealed y coordinate used in crystal generation 


JTS: 


Variable used to establish atom sites 


JTT: 


n n tt f! I! II 


KF: 


Final K in LOCAT (K) assigned to an atom 


KT: 


Unsealed z coordinate used to establish atom site 


LCUT( I) 


: Used to identify an ith atom which is not included in 
calculations 


LD: 


The highest numbered atom in the mobile layers 


LL: 


The highest numbered atom in the entire crystal 



LOCAT(K) : Dimensioned variable that remembers the numbers of the 





atoms within a radius ROEL of the primary at time zero 


LS: 


Variable associated with each of the nine lattice 
generator subr outines 


MCRO: 


One number higher than the order of the fit between the 
Born-Mayer and Morse potentials, always 4 in this simu- 
lation 


ND: 


Data output increment, in numbers of timesteps 


NEW: 


Parameter used to determine whether or not atom numbers 
have been stored in LOCAT (K) 


NPAGE: 


Page numbering variable 


NRUN: 


Parameter used to determine whether or not to read 




additional data cards 



57 



NS: 


Initial print statement timestep number 


NT: 


Timestep number 


NTT: 


Timestep number limit before shutdown 


NVAC: 


An atom number used to establish point defects or used as 



a reference point for interstitial placement 



PAG: 


Parameter for bullet force function correction 


PBMAS: 


Primary mass in kilograms 



PEXA, PEXB: Input Born-Mayer potential function parameters for 





the bullet “target interaction 


PFPTC: 


Primary force function evaluated at ROE 


PFXA: 


Primary force function parameter 



PKE(I): Kinetic energy of the i th atom 



PLANE : 


Alpha-numeric array for lattice orientation 


POT: 


Potential energy between two atoms 



PPE(I) : Potential energy of the i^ atom 

PPTC: Primary potential function evaluated at ROE 

PTE(I): Total energy of the i th atom (potential + kinetic) 



PTMAS: 


Target mass in kilograms 


RE: 


Input Morse potential parameter 


RO: 


Spacing constant in FGG(llO) lattice generation subroutine 


ROE: 


Nearest neighbor distance 


R0E2: 


ROE squared 


ROEA: 


Maximum cut off for Born-Mayer potential 


ROEB: 


Minimum cut off for Morse potential 


ROEG: 


Maximum cut off for Morse potential 


R0EG2 : 


ROEG squared 


ROEL: 


Radius inside of which local potential energy is i 'Und 
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R0EL2 : 


ROEL squared 


ROEM: 


ROE-DTI, region in which modification of repulsive force 




must be made 



RX(I), RY(I), RZ(I): x,y,z coordinates of an ith atom at any time 
RXI(I), RYI(I), RZI(I): x,y,z coordinates of an i th atom’s initial 
position 

RXK(I), RYK(I), RZK(I): x,y,z coordinates of temporary position of 





an i^ atom during force cycle 


SAVE: 


POT 



sex, SCY, SeZ: x,y,z coordinate scale factors 



SSCZ: 


A z scale factor used for the FCC(lll) lattice generator 



subr out ine 



START : 


An optional timing variable, not used in this simulation 


SUM; 


Variable in cubic fit subroutine 


TARGET: 


Alpha-numeric array for target material 


TSAVE : 


Bullet mass/( target mass + bullet mass); distributes 



potential energy between target and bullet 



TE: 


Total energy of all crystal atoms (kinetic + potential) 


TEMP: 


Temperature of lattice in degrees Kelvin. Not used in 



this simulation 



TFAC: 


A time factor ratio used to determine DT by maximum force 



method 



TFACB: 


TFAC for the bullet 


THERM: 


Thermal energy of atom. Not used in this simulation 


TIME: 


Elapsed problem time in seconds 


TLPE: 


Total local potential energy of atoms within a radius ROEL 


TLPE0: 


TLPE at time zero 
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TMAS: 



Target atom mass in amu 



TPKE: 


Total kinetic energy of all crystal atoms 


TPOT: 


Total potential energy of all crystal atoms 


VSS: 


Storage variable for velocity components 



VX(I), VY(I), VZ(I): x,y,2 components of ith atoms velocity 
X, Y, Z: Unsealed coordinates used in crystal generation 
YLAX(I): Relaxation in -y direction of layer in L.U. 

ZP : Floating point form of JTT 
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FIGURE 6. 

MOTION IN -Y DIRECTION (TOWARD SURFACE) 



j ^ ^ 




6A 



T I WESTER NUMBER 



FIGURE 7. 



BCC (100) ORIENTATION 
HARD SPHERE MODEL 



SOLID line: "y" PLANE 

DASHED line: " Y -j." PLANE 

INTERSTITIALS IN "y" PLANE: 

INT A IN <0I0> OPEN CHANNEL 
INT B HIDDEN FROM SURFACE 



2 LU 










\ 



BINDING ENERGY 




66 



V(R) 



F(R) 



FIGURE 9. 

POTENTIAL FUNCTION 
EROSION 



ACTUAL POTENTIAL 



ERODED POTENTIAL 
= V'(R) 




ROEA 



SEPARATION 



FORCE FUNCTION 
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c 

THIS PROGRAM GENERATES VARIOUS TYPES AND ORIENTATIONS OF 
CRYSTAL LATTICES, AND INJECTS A VACANCY, INTERSTITIAL, OR 
REPLACEMENT IMPURITY AT A DESIRED LOCATION. IT THEN, BY USE 
OF ATOMIC POTENTIAL PARAMETERS AND NEWTONIAN MECHANICS, CAL- 
CULATES THE DYNAMIC RESULTS OF THE SYSTEM; OUTPUTING POS- 
ITION, VELOCITY, AND ENERGY VALUES FOR EACH ATOM IN THE 
CRYSTAL. 

C 

DIMENSIONING OF VARIABLES NOT NEEDED IN COMMON 

DIMENSION VX(IOCO) , VY( lOOC ) ,VZ( lOOn) , PKE( lOCO) 
DIMENSION OX ( 1000 ) , DY( 1000 ) , DZ ( 1000 ) , PTE( 1000) 
DIMENSION RXK( ICOO) , RYK ( 1000 ), RZK ( 1000 ) 

C 

COMMON LABELING OF VARIABLES REQUIRED IN OTHER SUBROUTINES 
COMMON/ CCMl/RX ( 1000 ) , RYdOOO ) ,RZ( 1000 ) , LCUT( 1000 ) , 

ILL ,LD, ITYPE ,NVAC 

COMMON/ C0M2/ IHK 20) , I H2 ( 8 ) , I HS ( 10 ) , I HB ( 6 ) , I HT ( 6 ) , 
ITARGET (4) , TMA5 , BULLET(4) , BMAS , PL AN E , T EMP , T HERM 
C0MM0N/C0M3/RXI ( lOOC) , RYI ( 1000 ) ,RZI ( 1000) ,CVR,EVR, 
1NT,TIME, DT , DTI , ILAY 

C0MM0N/C0M4/IX ,I Y, I Z , SC X , SCY , SC Z , I DEE P , D1 X , D1 Y , D1 Z 
COMMON/COM5/ROE,ROE2,ROEM, E XA , E XB , P EXA , P E XB , F XA , P F XA , 
1IQ,TSAVE,BSAVE 

C0MM0N/C0M6/FX( 1000 ) ,FY( iCOO ) , FZ( 1000) , PAC,PFPTC, FM 
COMMON/C 0M7/PPTC,T POT, PPE( 1000 ) , TLPE,R0EL , R0EL2, NEW 
C0MM0N/C0M8/R0EA, ROFB, ROEC , R0EC2 , CPO , CP 1 , C P2 , CP3 , 
1CF0,CF1 ,CF2, CGD1,CGD2, CGBl ,CGB2,CGF1 ,CGF2 
COMMON/COMA/ A(4,5),MCR0 
C 

READ STATEMENT FORMATS 
9010 FORMAT (20A4) 

9020 FORMAT { 8A4,3F8. 5,2F5.2 ) 

9030 F0RMAT(4A4, 3F8. 5, 6A4,F6.2) 

9040 F0RMAT(F6.2,F5o3,I5,6I4,3F5.2, 12) 

90 50 FORMAT! U!A4,A4.4I3tF8.4,I4) 

r 

WRITE STATEMENT FORMATS 
9610 FORMAT! 1 HI) 

9620 FORMAT !47X, 'SUMMARY OF ATOMS ' / / , 3 5X , 8A4 , ' , NT =‘14,//, 
13C ATOM POSITION BIND ENERGY »),//) 

96 30 FORMAT! 3! I 5 , 3F6 . 2 , F 8. 4 , 8X ) ) 

9640 FORMAT! /4X,F1C. 3, 25H EV, TOTAL KINETIC ENE RGY, , F 10 . 3 , 
127H EV, TOTAL POTENTIAL ENERGY , FlC . 3 , 13H EV , RE DUCT I ON , 
1//,20X,F1U.3,50HEV, LOCAL POTENTIAL ENERGY, IN VOLUME 
lOF RADIUS = ,F5.2, /30X, 16HCHANGE IN TLPE =, FIG. 3) 

9650 FORMAT! 105X,4HPAGE , 13, /,1H1) 

9660 FORMAT!/ • ATOM DX DY DZ 

IVX VY VZ KE PE TE'/) 

9670 FORMAT! 118, 3F10c3,3F10. 1,3F10. 4 ) 

C 

INITIALIZING 

START = 0.01=4=ITI ME ! XX) 

DO 2 1=1,1000 
LCUT!I)=C 
RXK! I )=0.0 
RYK! I )=0.0 
RZK! I ) =0.0 
VX! I )=0.0 
VY! I )=0 .U 
VZ! I )=0.0 
PKE! I )=0.0 
PPE( I )=0.0 
PTE ! I ) =0.0 
RX! I )=G,0 
RY!I)=C.O 
RZ! I )=0.0 
FXI I )=0.0 
FY! I )=0.0 
FZ! I )=0.0 
RXl ! I ) =0.0 
RYI ! I )=0.0 
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2 RZK I )=0.0 
ISHUT=1 
NRUN=0 
C 

INPUT DATA 

READ ( 5,9010) 
READ ( 5,9C2D) 

READ ( 5,9030) 
READ ( 5,9G30) 
READ ( 5,9C40) 
1D1X,D1Y,D1Z, IQ 
READ ( 5,9050) 
C 

CONSTANTS AND SCALING 
R0E2=3.0 
ROE=SQRT(ROE2) 
ROEM = ROE-DTI 
ROEL2=ROEL^ROEL 



IHl 

IH2 ,DCON, ALPHA, RE ,ROEC, ROE L 
BULLET, BMAS,PEXA, PE XB, I HB, THERM 
TARGET,TMAS, EXA,EXB, 1HT,TEMP 
EVR,DTI ,NTT,NS,ND, IP, I DE E P , I TYPE , NV AC 

IHS,PLANE,LS, I X , I Y , I Z , CVR , MCRO 

FACTORS 



CVe=l*60E-19 
• EV=EVR^‘1.0E+3 
CVM=i.672E-27 
FM=1«0E-10 
FM2=FM*FM 
CVD=CVR*1.0E-10 
CVED=CVE/CVD 
PTMAS=THASYCVM 
PBMAS=BMAS-CVM 
HTMAS=0.5*PTMAS/CVE 
HBMAS = U« 5^'^PBMAS/CVE 
TSAVE=BMAS/( BMAS+TMAS) 
BSAVE=THAS/ ( BMAS+TMAS ) 



C 

REPULSIVE POTENTIAL PARAMETERS 
FXA=ALOG(-EXBACVED ) +EXA 
PFXA=ALOG( -PEXBYCVED) +PEXA 
PPTC = PX.P ( PFXA+PFXR*-!<RnE ) 

PAC-ALOG(CVED) +PEXA 
PFPTC=EXP( PAC+PEXB^ROE ) 

C 

ATTRACTIVE POTENTIAL PARAMETERS 

CGD1 = /.L0G( DCON ) +2. O^i'AL PHA*RE 

CGD2 = AL0G( 2.C*DC0N ) +ALPHA=i‘RE 

CGB1=-2.0*ALPHA«CVR 

CGB2=-ALPHA-CVR 

CGF1=ALOG(-CGB1«CVEO)+CG01 

CGF2=AL0G(-CGB2*CVED)+CGD2 

C 

CUTOFF DISTANCES FOR ATTRACTIVE AND REPULSIVE POTENTIALS 
ROEA = ] « 50/CVR 
R0EB=2,0/CVR 
R0EC2=R0ECX=R0EC 






C 

PARAMETERS FOR CALCULATION OF THE BEST CUBIC FIT IN THE GAP 
BETWEEN MAXIMUM DISTANCE CUTOFF OF THE REPULSIVE POTENTIAL 
(ROEA), AND MINIMUM DISTANCE CUTOFF OF THE ATTRACTIVE POTEN- 
TIAL (ROEB). SUBROUTINE CROSYM ACTUALLY PERFORMS THIS CURVE 
FITTING. 

A( 1,1)=1.0 

A( 1,2) = R0EA 

A(1,3)=R0EA>^R0EA 

A( 1,4) =R0EAX==^3 

A( 1,5) = EXP( EXA + EXB=<^ROEA) 

A(2,l)=1.0 

A(2,2)=ROEB 

A(2,3)=R0EB=^=R0EB 

A(2,4)=R0EB*-:=3 

A(2,5) = EXP(CGD1+CGB1*R0EB)-EXP(CGD2+CGB2=^'R0EB) 

A( 3,1 )=0.0 

A(3,2)=-1.0 

A(3,3)=-2.0=^RDEA 

A( 3,4) =-3.0*R0EA’;^R0EA 

A(3,5) =EXP( FXA+EXB^ROEA)/CVED 
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A(4, 1 )=0.0 
A(4,2)=-1.C 
A( 4,3) =-2.0=!=ROEB 
A(4,4)=-3.CvR0EB’!'R0EB 

A(4, 5) = ( EXP(CGE1+CGB1*R0EB)-EXP (CGF2+CGB2^'ROEB) ) /CVED 
CALL CROSYM 
CPO=A( 1 ,5) 

CP1=A(2,5) 

CP2=A(3,5) 

CP3=A(4, 5) 

CFO=-CPl=!=CVED 
CF1=-2.0»^=CP2=!^C VED 
CF2=-3.0=^'CP3YCVED 
C 

SELECTION OF THE DESIRED CRYSTAL STRUCTURE AND ORIENTATION, 
100, 110, AND 111 PLANES OF FACE-CENTERED, BODY-CENTERED, 
AND DIAMOND STRUCTURES ARE ALLOWED. I L AY AND IDEEP ARE VAR 
lABLES ESTABLISHING THE NUMBER OF MOBILE LAYERS IN THE 
CRYSTAL. P.XI(I) AMD RXK ( I ) ARE VARIABLES SAVING THE ORIGIN 
AL X-POSITION OF THE I’TH ATOM. Y AND Z POSITIONS ARE 
ANALOGOUS. 

GOTO ( 11 , 12, 13, 14,15 , 16, 17, 18, 19) ,LS 

11 CALL LlOO 
GO TO 30 

12 CALL 1.110 
GO TO 30 

13 CALL Llll 
GO TO 30 

14 CALL BlOO 
GO TO 3G 

15 CALL BUG 
GO TO 30 

16 CALL Bill 
GO TO 30 

17 CALL DICO 
CQ TO ^0 

18 CALL DllG 
GO TO 30 

19 CALL Dill 

30 ILAY=IOEEP 

IF (IDEEP) 35,35,40 

35 LD=LL 
ILAY=IY 

40 DO 45 I=1,LL 
RXK( I )=RX( I ) 

RYK( I )=RY( I ) 

RZK( I ) =RZ( I ) 

RXI ( I )=RX( I ) 

RYI ( I )=RY( I ) 

45 RZI ( I ) =RZ( I ) 

C 

THIS SECTION ALLOWS (JNE TO REPEAT A RUN OF THE PROGRAM WITH 
DIFFERENT DATA WITHOUT REPEATING INITIALIZATION, POTENTIAL 
PARAMETER CALCULATIONS AND CRYSTAL LATTICE BUILDING. SUB- 
ROUTINE PLACE USES LCUT(I) AND NVAC TO CREATE VACANCIES, 
INTERSTITIALS, AND REPLACEMENT IMPURITIES AT DESIRED LOCA- 
TIONS IN THE LATTICE, 

IF(NRUN.EQ.O) GO TO 60 

50 READ ( 5,9G40) E VR ,DT I , NTT , NS , NO , I P , I DEE P , I TYPE, NV AC 
1D1X,D1Y,D1Z, IQ 
IF(DTI.EO.O) GO TO 9999 
DO 55 1=1, LL 
LCUT( I )=G 
RX( I ) =RXI ( I ) 

RY( I )=RYI ( I ) 

RZI I ) = RZI ( I ) 

RXK( I ) =RXI ( I ) 

RYK( I )=RYI ( I ) 

55 RZK( I ) =RZI ( I ) 

60 NRUN=1 

CALL PLACE 
RXI(1)=RX(1) 
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I 
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RYI (1)=RY( 1) 

RZI ( 1) = RZ( 1 ) 

RXK(1 )=RX( 1 ) 

RYK( 1)=RY( 1) 

RZK(1 ) = RZ( 1 ) 

DO 65 1=1, LL 
VX( I )=C.O 
VY( I )=G.O 
VZ( I )=0.0 
PPE( I ) =0.0 
PKE( I )=C,0 
PTE( I )=0.0 
TP0T=0.0 
NEW=0 
C 

THE ENERGY SUBROUTINE CALCULATES THE POTENTIAL ENERGY OF 
EACH ATOM IN THE LATTICE. SUBROUTINE LOCAL SUMS UP THIS 
ENERGY FOR ALL ATOMS WITHIN A SPECIFIC RADIUS OF THE POINT 
DEFECT . 

CALL ENERGY 
CALL LOCAL 
BINO=~TPOT 
TPKE=0.0 
TLPEO=TLPE 
DLPE=TLPE-TLPEO 
TE=TPOT+BIND 
C 

THIS SECTION PRINTS OUT X, Y, AND Z COORDINATES, IN LATTICE 
UNITS, AND BINDING ENERGIES OF EACH ATOM IN THE CRYSTAL AT 
TIME ZFRO. 

TIME=0.0 

NT=0 

WRITE ( 6,9610) 

WRITE ( 6,9620) IH2,NT 
DO 70 1=1, LL, 3 
K=I + 1 
J = I ■’•' 2 . 

70 WRITE ( 6,9630) I , RX ( I ) , RY ( I ) , RZ ( I ) , PP E U ) , K , RX ( R ) , 

1RY(K),RZ(K),PPE(K),J,RX(J),RY(J),RZ(J),PPE(J) 

WRITE ( 6,96A0) T PKE , T POT , TE , TL P E , ROEL , D L PE 
NPAGE=1 
NPAGE=NPAGE+1 
WRITE ( 6,9650) NPAGE 
C 

THIS IS THE MAIN BODY OF THE PROGRAM. BY USE OF THE AVERAGE 
FORCE METHOD, EXPLAINED IN DETAIL IN APPENDIX C, IT DOES ALL 
THE DYNAMICS FOR EACH INDIVIDUAL ATOM. SUBROUTINE STEP 
CALCULATES ALL MUTUAL FORCES AMONG THE ATOMS. BASED ON THE 
FORCES, THIS SECTION THEN CALCULATES TEMPORARY POSITIONS FOR 
THE PRIMARY, AND ALL OTHER ATOMS; RECALCULATES FORCES IN 
STEP; AND THEN RECALCULATES FINAL POSITIONS FOR THE PRIMARY 
AND ALL OTHER ATOMS, BASED ON THE AVERAGE OF THESE TWO 
FORCES. THIS SECTION ALSO INCLUDES ALL KINETIC ENERGY CAL- 
CULATIONS, BASED ON THE VELOCITIES INVOLVED; AND FINALLY 
CALCULATES A NEW TIMESTEP DURATION FOR USE IN THE NEXT TIME- 
STEP, BASED ON EITHER A MAXIMUM ALLOWED FORCE, OR MAXIMUM 
ALLOWED ENERGY. (SEE APP. C) VELOCITIES ARE HALVED AT THE 
END OF EACH TIMESTEP AS A METHOD OF DAMPING. 

95 TFAC = 2.0=:=PTMAS=!‘DTI*CV0 
TFACB=2.0={'PBMAS=!=DTI*CVD 
DT=1.0E-14 
100 DT09=DT/CVD 

HDT0D=C.5=i'DT0D 
DTOM=OT/PTMAS 
HDT0H=0. 5*DT0M 
DTOHB=DT/PBMAS 
HDT0MB=0.5*DT0MB 
200 CALL STEP 

IF(LCUT( D.GT.O) GO TO 240 
1=1 

RXK( I )=RX( I ) 

RYK( I )=RY( I ) 
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RZK( I ) =RZ( I ) 

RX( 1 )=RX( I )+DTOD«( HDTOMB*FX( I )+VX( I ) ) 

RY( I)=RY(I )+DTOD=!=( HDT0M8*FY( I)+VY( I) ) 

RZ( I ) = RZ( I )+OTQD’!M HDTOMB*FZ( I )■» VZ( I ) ) 

240 DO 245 1=2, LO 

IF(LCUT( n .GT.O )G0 TO 245 
RXK( I ) =RX( I ) 

RYK( I )=RY( I ) 

RZK( I )=RZ( 1 ) 

RX( I )=RX(I )+OTOD*( HDTOM*FX( I )+VX( I ) ) 

RY( I ) = RY( I )+DTOD-'!^( HDTOM-FY( I ) + VY( I ) ) 

RZ( I )=RZ( I ) +DTOD*( HOTOM*FZ( I )+VZ( I) ) 

245 CONTINUE 
CALL STEP 
EMAX=0,0 
FMAX=0. 0 
TIME=TIME+DT 
NT=NT+1 

IF(LCUT( 1) .GT.O) GO TO 265 
1 = 1 

VSS=VX( I ) 

VX( I )=VSS+HDTOMB=!'FX ( I ) 

RX( I )=RXK( I )+( VX( I )+VSS )=.’'HDTOD 
VSS=VY( 1 ) 

VY( I )=VSS+HDTOMB=^'FY(I ) 

RY( I)=RYK(I )+(VY(I )+VSS)*HDTOD 
VSS=VZ( I ) 

VZ( I )=VSS+HDTOMB*FZ( I ) 

RZ( I ) = RZK( I ) + ( VZ ( I )+VSS )''.''HDTOD 

PKE( I )=VX( I )*VX( I )+VY( I )*VY( I ) + VZ ( I ) « VZ ( I ) 

EMAX=PKE( I ) 

FMAX=FX( I )=;=FX( I )+FY(I )^FY( I ) + FZ ( I )=i^FZ ( I ) 
F0RC1=FMAX 
260 FX(I)=0«0 
FY( I )=e.O 
FZ( I » = o.n 
265 DO 280 1=2, LD 

IF(LCUT( I)oGT.0)G0 TO 280 
VSS=VX( I ) 

VX( I)=VSS+HOTOM*FX( I ) 

RX( I ) = RXK( I ) + ( VX( I )+VSS)^=HDTOD 
VSS=VY( I ) 

VY( I ) =VSS+HDTOM'^FY( I ) 

RY( I ) = RYK( I ) + ( VY( I )+VSS )=^HDTOD 
VSS=VZ( I ) 

VZ( I ) =VSS+HDTO/-K'FZ( I ) 

RZ( I )=PZK( I )+( VZ ( I ) +VSS )-HOTOD 
PKE ( I ) =VX( I )*VX( I )+VY{ I )>"VY( I ) +VZ ( I )--!-'VZ ( I ) 
2 75 F2 = FX( I )^FX( I )+FY( I ) *F Y ( I ) +FZ ( I ) =<=F Z ( I ) 

FX( I )=0.0 
FY(I)=C.O 
FZ( I )=C .0 

IF(F2<.GToFMAX) FMAX = F2 
IF(PKE( I I.GT.EHAX) EMAX=PKE( I ) 

280 CONTINUE 

I F( EMAX. EQ.O.G ) GO TO 285 
IF(FMAX,EQ.U.O) GO TO 285 
GO TO 287 
285 DT=1.0E-14 
GO TO 500 

287 IF(EMAX.EQ.PKE( 1 ) ) GO TO 290 
IF(FHAX.EQ. FORCl) GO TO 290 
EMAX=SQRT( EMAX ) 

FMAX=SQRT( FHAX ) 

DT=SQRT (TFAC/FMAX) 

FTERN=FMAX^DT/ ( 2.0*PTMAS) 

GO TO 295 

290 EMAX=SQRT(PKE( 1) ) 

FMAX=SORT( FORCl ) 

IF ( FORCl .EO.O.C ) GO TO 285 
DT=SQRT(TFACB/FMAX) 

FTERK=FMAX*DT/ (2.0^PBMAS) 
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I 



A 



295 IF(EMAX.LE.FTERM) GO TO 300 
FDT I=DT I*CVD 
DT=FDTI /EMAX 

300 IF ( ISHUT.EQ.-l ) GO TO AOO 
310 IF(NS-NT) 400f 400,320 
320 DO 325 1=1, LL 
VX( I )=(). 5=!'VX( I ) 

VY( 1 )=0 ,5-VY( I ) 

325 VZ( I )=G.5-VZ(I ) 

GO TO 100 
340 DO 350 I=1,LL 
RX( I )=RXK( I ) 

RY( 1 )=RYK( I ) 

350 RZ( I )==RZK( I ) 

NTT=NT 

C 

THE PRINT SUBROUTINE PLACES A HEADING OF PERTINENT INFORMA- 
TION AT THE TOP OF EACH TIMESTEP PRINTOUT. 

400 CALL PRINT 
C 

POTENTIAL ENERGY AND LOCAL POTENTIAL ENERGY FOR EACH ATOM 
ARE CALCULATED BASED ON THE NEW POSITIONS. SUMMATIONS OF 
TOTAL POTENTIAL AND KINETIC ENERGY FOR THE LATTICE ARE PER- 
FORMED. DX, DY, AND DZ KEEP TRACK OF MOTION RELATIVE TO THE 
INITIAL POSITION AT TIME ZERO FOR EACH ATOM. 

410 TP0T=0.0 

DO 450 I =1 , LL 
PPE( I )=0.0 
450 PTE( n =C.O 
CALL ENERGY 
CALL LOCAL 
DLPE=TLPE-TLPEO 
PKE(1)=HBMAS^PKE(1) 

TPKE=PKE(1) 

PTE(1)=PKE(1)+PPE( 1) 

DO 620 1=2, LL 

PKF f T ) =MTMAS-PKF ( T 1 

TPKE=TPKE+PKE( I ) 

620 PTE( I)=PKE( I)+PPE( I) 

TE=TPOT+BIND 
WRITE ( 6,966D) 

700 DO 750 1=1, LD 

0X( I )=RX( I )-RX Id) 

DY( I)=RY( I )-RYI ( I ) 

DZ( I) = RZ( I )-RZ Id) 

C 

THIS SECTION PRINTS THE RELATIVE MOTION, VELOCITY, AND 
ENERGY OF EACH ATOM, FOR EVERY TIMESTEP SO DESIGNATED: IE, 

EVERY ND'TH TIMESTEP, BEGINNING WITH #NS AND ENDING WITH 
#NTT, 

720 WRITE ( 6,9670) I , DX ( I ) , DY { I ) , D Z d ) , VX ( I ) , VY ( I ) , 

1VZ( I ) ,PKE( I ) , PPE( I ) , PTE( I ) 

750 CONTINUE 

WRITE ( 6,9640) T PKE , T POT , T E , T LP E , ROEL , D L PE 
WRITE ( 6,9650) NPAGE 
NPAGE=NPAGE+1 
IF(NT-NTT) 760,950,950 
760 DO 780 1=1, LL 
VX( I )=0.5^VX( I ) 

VY( I )=(). 5*VY( I ) 

VZ( I )=0.5«VZ( I ) 

780 CONTINUE 
790 NS=NS+ND 
GO TO 100 
950 CONTINUE 
C 

THIS SECTION PRINTS OUT X, Y, AND Z COORDINATES AND BINDING 
ENERGIES OF EACH ATOM IN THE CRYSTAL AT THE END OF THE 
PROGRAM. 

955 WRITE ( 6,9620) IH2,NT 
DO 965 1=1 , LL, 3 
K=I + 1 
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J=I+2 

965 WRITE ( 6,9630) I , RX ( I ) , RY ( I ) , RZ ( I) , PP E ( 1 ) , K , RX( K ) , 
1RY(K),RZ(K),PPE(K),J,RX(J),RY(J),RZ(J),PPE(J) 

WRITE ( 6,9640) TPKE , TPOT , TE , TLPE , RCEL , DLPE 
WRITE ( 6,9650 NPAGE 
1000 IF(ISHUT) 9999,50,50 
9999 STOP 
END 



SUBROUTINE CROSYM 
C 

SOLVES R SIMULTANEOUS EQUATIONS BY THE METHOD OF GROUT 
THIS SUBROUTINE FITS THE BEST CUBIC BETWEEN THE REPULSIVE 
AND ATTRACTIVE PARTS OF THE POTENTIAL. 

C 

COMMGN/COMA/ A(4,5),MCR0 
M=MCRO 
N = M+1 
1 1 = 1 

100 13=11 

SUM=ABS(A( 11,11)) 

DO 120 I=I1,M 

IF(SUM-ABS (A( I , ID ) ) 110,120,120 
110 13=1 

SUM=ABS ( A( I , II ) ) 

120 CONTINUE 

IF(13-I1) 130,150,130 
130 DO 140 J=1,N 
SUM=-A( II, J) 

A( II, J)=A( 13, J ) 

140 A( 13, J ) = SUM 
150 13=11+1 

DO 160 I=I3,M 

160 A( I , II ) = A( I , II )/A( I 1, I 1 ) 

170 J2=T1-1 

I T 14- 1 

IF(J2) '180,200, 180 
180 DO 190 J=I3,N 
DO 190 1=1, J2 

190 A( II, J)=A( I1,J )-A( II, I )«A( I, J) 

IF(Il-M) 200,220,200 
200 J2=I1 
11 = 11+1 
DO 210 I=I1,M 
DO 210 J=1,J2 

210 A( I , 11 )=A( I , II ) - A( I ,J)=!'A( J, I 1) 

IF(Il-M) 100,170,100 
2 20 DO 240 1 = 1, M 
J2=M-I 
I3=J2+1 

A( I3,N)=A( I3,N)/A( 13,13) 

IF(J2) 230,250,230 
230 DO 240 J=1,J2 

240 A(J,N)=A(J,N)-A(I3,N)#A(J,I3) 

250 RETURN 
END 



SUBROUTINE LlOO 
C 

THIS IS A LATTICE GENERATOR FOR THE FCC (100) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, Z FOLLOWED BY Y, 
FOLLOWED BY X. 

IT CONTAINS A NONSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

COMMON/COMl/RX ( 1000 ) , RY( 1000 ) , RZ( 1000 ) , LCUT ( lOGC ) , 
ILL, LD, ITYPE ,NVAC 

C0MM0N/C0M4/1 X, lY, I Z , SCX , SCY , SC Z , IDEEP,D1X,D1Y,D1Z 
DIMENSION YLAX(20) 

DATA YLAX/2C*0.0/ 
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YLAX( 1 )=-0.097 

YLAX(2) =-0.024 

SCX=1.0 

SCY=1.0 

SCZ=1.0 

M=2 

JT=0 

Y=-SCY 

DO 60 J=1,IY 
Y=Y+SCY 
KT=0 
Z=-SCZ 

DO 59 K=l, IZ 
Z=Z+SCZ 
IT=0 
X=-SCX 

DO 58 1=1, IX 

X=X+SCX 

ITT=IT+JT+KT 

IF( ITT-( ITT/2)=!‘2) 57, 30,57 
30 RX(M>=X 

rY(M)=Y+YLAX(J ) 

RZ(M)=Z 

M=M+1 

57 CO.'^TINUE 
IT =IT+1 

58 CONTINUE 
KT=KTH 1 

59 CONTINUE 

JT = JT + 1 

IF( IDEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 

100 RETURN 
110 LD=M-1 
GO TO 60 
END 



SUBROUTINE LllO 
C 

THIS IS A LATTICE GENERATOR FOR THE FCC (110) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, Z FOLLOWED BY Y, 
FOLLOWED BY X. 

IT CONTAINS A NONSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

C0MM0N/C0M1/RX( 1000 ), RY( 1000 ) ,RZ( 1000 ) ,LCUT( lOOC ) , 
1LL,LD, ITYPE,NVAC 

C0MM0N/C0M4/IX, lY, I Z , SC X , SCY , SC Z , I DEEP , D1 X , D1 Y , D1 Z 
DIMENSION YLAX(20) 

DATA YLAX/20^0.0/ 

YLAX( 1 ) =-0.07 
YLAX(2 ) =-0.02 
R0=1.0/SQRT( 2. 0) 

SCX=R0 

SCY=RO 

SCZ=1.0 

M=2 

JT=0 

Y=-SCY 

DO 60 J=1,IY 
Y=Y+SCY 
KT=0 
Z=-SCZ 

DO 59 K = l, IZ 
Z=Z+SCZ 



DO 58 1 =1, IX 

Y= y4- y 

IF( lT-( IT/2)*2 > 21,11,21 
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11 IF( JT- ( JT/2 )*2 ) 57,12,57 

12 IF(KT-(KT/2 )=>=2 ) 57,30,57 

21 IF( JT-( JT/2 )=^=2 ) 22,57,22 

22 IF( KT-( KT/2 )=i‘2 ) 30,57,30 
30 RX(M)=X 

RY(,M) =Y + YLAX( J) 

RZ(M)=Z 

M=M+1 

57 CONTINUE 
IT =IT+1 

58 CONTINUE 
KT=KT+1 

59 CONTINUE 

JT = JT + 1 

IF(IOEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 

100 RETURN 
110 LD=M-1 
GO TO 60 
END 



SUBROUTINE till 
C 

THIS IS A LATTICE GENERATOR FOR THE FCC (111) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, Z FOLLOWED BY Y, 
FOLLOWED BY X. 

IT CONTAINS A. NONSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

COMMON/ CCM1/RX( 1000 ) , RY( IGOD ) , R Z (1000 ) , LCUT ( 1000) , 
ILLtLD, ITYPE,NVAC 

C0MM0N/C0M4/IX, lY, I Z , SCX , SCY , S C Z , I DEEP , D 1 X , D1 Y , D 1 Z 
DIMENSION YLAX(2C) 

DAT A Y LAX/ ? , 0 / 

YLAX( 1 ) =-0.0A 
YLAX(2 )=-C.01 
SCX=1.0/SQRT(2.0) 

SCY=2«(;/SQRT(3.0) 

SCZ=SQRT( 1. 5) 

SSCZ=SCZ/3cO 

H=2 

JT=0 

Y.— — 

DO 60‘ J=l, lY 

Y=Y+SCY 

JTS=JT+JT/3 

Z=-SCZ 

KT=0 

DO 59 K=l, I Z 
Z=Z+SCZ 
IT=0 
X=-SCX 

00 58 1 = 1, IX 

X=X+SCX 

IN=IT+JTS+KT 

IF( IN-( IN/2 )=f-2 ) 57,30,57 
30 RX(M)=X 

RY(M)=Y + YLAX( J ) 

IF( JT-3*( JT/3) ) 41,45,41 

41 JTT=JT 

42 JTT=JTT-3 
IF(JTT) 43,45,42 

43 JTT=JTT+3 
ZP=JTT 

RZ(M)=Z+ZP«SSCZ 
GO TO 50 
45 RZ(M)=Z 
50 H=M+1 
57 CONTINUE 
IT=IT+1 
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58 CONTINUE 
KT = KH-1 

59 CONTINUE 
JT=JT+1 

IF(IDEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 

100 RETURN 
110 LD=M-1 
GO TO 60 
END 



SUBROUTINE BlOO 
C 

THIS IS A LATTICE GENERATOR THE THE BCC (100) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, X FOLLOWED BY Z, 
FOLLOWED BY Y. 

IT CONTAINS A NONSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

C0MH0N/C0M1/RX( lOCO ) ,RY( 1000 ) ,RZ( 1000) ,LCUT( lOOC) , 
1LL,LD, ITYPE,NVAC 

COMMON/ COMA/ I X, I Y. I Z, SCX, SCY , SCZ , I DEEP , Dl X ,D1Y, D1 Z 
DIMENSION YLAX (20) 

DATA YLAX/20=^P.0/ 

YLAX( l)=-0c2G 

YLAX(2)=-0.03 

SCX=1.0 

SCY=1,0 

SCZ=1 .0 

M = 2 

JT=0 

Y=-SCY 

DO 60 J=l, I Y 
Y=Y+SCY 

t/ T-ri 

Z=-SCZ 

DO 59 K=1,IZ 

z--=z+scz 

IT=0 

X=-SCX 

DO 58 1=1, IX 
y = y+ <^r y 

IF(IT-( IT/2)*2) 21, 11,21 

11 IF( JT-( JT/2 )=^2 ) 57,12,57 

12 IF(KT-(KT/2)-2) 57,30,57 

21 IF( JT-( JT/2 )-2) 22,57,22 

22 IF( KT-(KT/2)=:‘2 ) 30,57,30 
30 RX(M)=X 

RY(M)=Y+YLAX( J ) 

RZ(M)=Z 
M=M + 1 

57 IT=IT+1 

58 CONTINUE 
KT=KT+1 

59 CONTINUE 

JT = JT + 1 

IF(IDEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 

100 RETURN 
110 LD=M-1 
GO TO 60 
END 
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SUBROUTINE BllO 
C 

THIS IS A LATTICE GENERATOR THE THE BCC (IIC) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, X FOLLOWED BY Z, 
FOLLOWED BY Y. 

IT CONTAINS A NONSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

COMMOM/CCMl/RX( 1000 ) , RY( 1000 ) , RZ( lOOC ) , LCUT( 1000 ) , 
ILL, LD, ITYPE,NVAC 

COMMON/ COMA/ IX, lY, I Z , SCX , SCY , SCZ , I DEEP , D1 X , D1 Y , D1 Z 
DIMENSION YLAX(2G) 

DATA YLAX/20*C.0/ 

YLAX( 1 ) =-0.10 
YLAX( 2 )=-0.01 
SCX = SQRT(2.D) 

SCY = SQRT(2.0) 

SCZ=1.C 

M=2 

JT=0 

Y=-SCY 

DO 60 J=l, lY 
Y=Y+SCY 
KT=0 
Z=-SCZ 

DO 59 K=l, IZ 
Z=Z+SCZ 
IT=0 
X=-SCX 

DO 58 1=1, IX 

X=X+SCX 

ITT=IT+JT+KT 

IF( ITT- ( ITT/2) *2) 57, 30,57 

30 RX(M)=X 

RY ( M) =v+vi AX ( J ) 

RZ(M)=Z 
M = M+1 

57 CONTINUE 
IT =IT+1 

58 CONTINUE 
KT=KT+1 

59 CONTINUE 

JT = JT + 1 

IF(IDEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 
RETURN 

no LD = M-1 
GO TO 60 
END 



SUBROUTINE Bill 
C 

THIS IS A LATTICE GENERATOR FOR THE BCC (111) ORIENTATION. 
THE CRYSTAL IS DEVELOPED IN THE ORDER, X FOLLOWED BY Z, 
FOLLOWED BY Y. 

IT CONTAINS A NGNSTANDARD USE OF THE SURFACE RELAXATION 
PARAMETER. 

C 

C0HM0N/CGM1/RX( ICOO ) ,RY( 1000 ) ,RZ( 1000 ) , LCUT ( 1000) , 
ILL, LD, ITYPE ,NVAC 

COMMON/ COMA/ IX, lY, I Z , SCX , SCY , SCZ , IDEEP,D1X,D1Y,D1Z 
DIMENSION Y LAX (20) 

DATA YLAX/2C’^0.0/ 

YLAX( 1 )=-0. 10 
YLAX(2) =-0.01 
SCX = SQRT(2.0) 
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SCY = SQRTd. 0/3.0) 

SCZ = SORT (6.0) 

M = 2 

JT=0 

Y=-SCY 

DO 60 J=l, I Y 
Y=Y+SCY 
KT = 0 

JJ = J-(JT/3)*3 
GO TO ( 5,10 ,15 ) , JJ 
5 Z =-SCZ 
GO TO 20 

10 Z =-SCZ+SCZ/3.0 
GO TO 20 

15 Z =-SCZ+(2.0>^SCZ/3.0) 

20 DO 59 K=1,IZ 
Z=Z+SCZ 
IT = 0 
X=-SCX 

DO 58 1=1, IX 
X=X+SCX 
ITT=IT+JJ+KT-1 
IF(ITT-( ITT/2) *2) 57,30,57 
30 RX(M)=X 

RY( M) =Y+YLAX( J ) 

RZ(M)=Z 
M = M+1 

57 IT=IT+1 

58 CONTINUE 
KT=KT+1 

59 CONTINUE 

JT = JT + 1 

IF(IDEEP-JT) 60,110,60 

60 CONTINUE 
LL=M-1 
RETURN 

110 LD=M-1 

GO TO 60 
END 



SUBROUTINE DlOO 

C0MM0N/C0M1/RX( 1000 ), RY ( 1000 ) ,RZ(1000 ) , LCUT(IOCO) , 
1LL,LD: ITYPE,NVAC 

COMMON/CCMA/IX, lY, IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z 

RETURN 

END 



SUBROUTINE DUO 

C0MM0N/C0M1/RX( 1000) ,RY( 1000) ,RZ( 1000) ,LCUT( 1000) , 
1LL,LD, I TYPE, NV AC 

C0MM0N/C0M4/I X, I Y, I Z , SCX , SCY , SC Z , I DE EP , D1 X , D1 Y , D1 Z 

RETURN 

END 



SUBROUTINE Dill 

C0MM0N/C0M1/RX( 1000 ) ,RY{ 1000 ) ,RZ( 1000) , LCUT( 1000) , 
1LL,LD, I TYPE, NV AC 

C0MM0N/C0N4/I X, I Y, I Z , SCX , SCY , SCZ , IDEEP,D1X,D1Y,D1Z 

RETURN 

END 



SUBROUTINE PLACE 
C 

THIS SUBROUTINE LOCATES A VACANCY, INTERSTITIAL, OR REPLACE- 
MENT IMPURITY IN THE LATTICE. 

C 

COMMON/COM1/RX( lOOu ) , RY( 1000 ) ,RZ(1C00 ) , LCUT( 1000) , 
ILLjLD, ITYPE,NVAC 
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C0MM0.\'/CGM4/IX , lY, IZ,SCX,SCY,SCZ,IDEEP,D1X,D1Y,D1Z 
GO TO { 10,20,30), ITYPE 
10 LCUT(NVAC) = 1 
LCUT(l) = 1 
RX( 1)=0.0 
RY( 1 )=0.0 
RZ( 1)=0.0 
GO TO AO 



20 


RX( 1 ) 


= 


RX ( NVAC) 


- DIX 




RY( 1) 


=r 


RY( NVAC) 


- DIY 




RZ( 1 ) 


= 


RZ( NVAC) 


+ DIZ 




GO TO 


40 






30 


LCUT(NVAC) = 1 






RX( 1) 


= 


RX ( NVAC) 






RY( 1) 


— 


RY(NVAC) 






RZ( 1) 


= 


RZ(NVAC) 




40 


RETURN 







END 



SUBROUTINE STEP 
C 

THIS SUBROUTINE DOES THE DYNAMICS FOR ONE TIMESTEP. 

THE FIRST HALF DOES THE DYNAMICS FOR ATOM «1; THE SECOND 
HALF FOR ALL OTHERS. 

C 

COMMON/ CCMl/RX ( 1000 ) ,RY( 1000 ) ,RZ( 1000) , LCUT( 1000) , 
ILL, LD, ITYPE, NV AC 

COMMON/CCM5/ROE,ROE2,ROEM,EXA,EXB, PEXA, PEXB,FXA,PFXA, 
1IQ,TSAVE, BSAVE 

C0MM0N/CGM6/FX ( lOCO ) ,FY( lOOU ) ,FZ( 1000 ) , PAC, PFPTC, FM 
C0MM0N/C0M8/R0EA, ROEB, ROEC,ROEC2,CPO,CP1 ,CP2,CP3, 
1CF0,CF1 , CF2 ,CGD1 ,CGD2,CGB1 ,CGB2,CGF1,CGF2 
IF( IQ.EQ.l) GO TO 200 
1 = 1 

!F{LCUT(I) ) 200, 105, 200 
105 I ^ + ’ 

" " DO i95 J=IP,LL 

IF(LCUT(J)) 195,110,195 
1 10 DRX = RX( J )-RX( I ) 

IF(DRX) 113,117,117 
113 IF(DRX+ROE) 195,195,120 
117 IF(DRX-ROE) 120,195,195 
120 DRY=RY( J)-RY(I ) 

IF(DP.Y) 123,127,127 
123 IF(DRY+ROE) 195,195,130 
127 IF(DRY-ROE) 130,195,195 
130 DRZ=RZ( J )--RZ( I ) 

IF(DRZ) 133,137,137 
133 IF(DRZ+ROE) 195,195,140 
137 IF(DRZ-ROE) 140,195,195 
140 DIST = DRX=^=DRX+DRY'-l=DRY+DRZ=:<DRZ 
IF(DIST-R0E2) 150,195,195 
150 DIST=SQRT( DIST ) 

160 IF(DIST-ROEM) 162,162,165 
162 FORCE=EXP( PFXA+PEXB*DIST) 

GO TO 180 
165 DFF=ROE-DIST 

IF( DFF-l.OE-10 ) 195,195,167 

167 FORCE=( EXP( PAC + P EX B*D I ST ) - PFPTC ) / DFF 
180 IF(FM-FORCE) 190,190,195 
190 FCD=FORCE/DIST 
FA = FOD=i^DRX 
FX( J)=FX( J )+FA 
FX( I ) = FX( I )-FA 
FA=FOD=<=DRY 
FY( J )=FY( J ) +FA 
FY( I )=FY( I )-FA 
FA=FOD*DRZ 
FZ( J ) = FZ( J )+FA 
FZ( I ) = FZ( I )-FA 
195 CONTINUE 



So 



c 



200 DO 300 I=IQ,LD 

IF(LCUT( I) ) 300,205,300 
205 IP=I+1 

DO 295 J=IP,LL 
IF(LCUT( J) ) 295,210,295 
210 DRX=RX( J)-RX( I ) 

IF(DRX) 213,217,217 
213 IF(DRX+ROEC) 295,295,220 
217 IF(DRX-ROEC) 220,295,295 
220 DRY=RY( J )-RY( I ) 

IF(DRY) 223,227,227 
223 IF(DRY+ROEC) 295,295,230 
227 IF(DRY-ROEC) 230,295,295 
230 DRZ=RZ( J )-RZ( I ) 

IF(DRZ) 233,237,237 
233 IF(DRZ+ROEC) 295,295,240 
237 IF(DRZ-ROEC) 240,295,295 
240 D1ST = DRX=^DRX+DRY^=DRY+DRZYDRZ 
IF(DIST-R0EC2) 250,295,295 
250 DIST=SORT(DIST ) 

IF( DIST-ROEA) 260,255,255 
255 IF(DIST-ROEB) 265,270,270 
260 FORCE=EXP(FXA+EXB=!'DIST ) 

GO TO 230 

2 65 FORCE = DI ST^i'CDI ST^^CF 2+C F 1 ) +CFO 
GO TO 280 

2 70 FORCE = EXP(CGF1+CG31YDI ST ) -EXP ( CGF2+CGB2=«'DI ST ) 
280 IF(ABS( FORCE), LE.FM) GO TO 295 
FOD = FORCE/DIST 
FA=F0D=SDRX 
FX( J )=FX( J )+FA 
FX( I )=FX( I )-FA 
FA=FOO*DRY 
FY( J )-FY(J )+FA 
FY( I )=FY( 1 )-FA 
PA=FQQ*nR7 
FZ( J ) = FZ( J )+FA 
FZ( I ) = FZ( I )-FA 
295 CONTINUE 
300 CONTINUE 
RETURN 
END 



SUBROUTINE ENERGY 
C 

THIS SUBROUTINE CALCULATES THE MUTUAL POTENTIAL ENERGIES. 
THE FIRST HALF DOES THE DYNAMICS FOR ATOM Ul; THE SECOND 
HALF FOR ALL OTHERS. 

C 

C0MM0N/C0M1/RX( lOCO) ,RY{ 1000 ) ,RZ( 1000) , LCUT( lOOC) , 
ILLjLD, ITYPE,NVAC 

COMMON/COM5/ROE,ROE2,ROEM,EXA,EXB, PEXA, PEXB,FXA,PFXA, 
1 IQjTSAVE, BSAVE 

COMMON/C GM7/PPTC ,T POT , PPE ( 1000 ) , TLPE , ROEL , R0EL2 , NEW 
C0MM0N/C0M8/R0EA,R0EB,R0EC,R0EC2,CP0,CP1 ,CP2 ,CP3, 
1CF0,CF1 ,CF2 ,CGD1 ,CGD2, CGBl , CGB2 , C GFl , CGF2 
IF(IQ.EQ.l) GO TO 20(J 
1 = 1 

IF(LCUT( I ) ) 600,505,600 
505 1P=I+1 

DO 595 J=IP,LL 
IF(LCUT(J)) 595,510,595 
510 DRX=RX( J )-RX( I ) 

IF(DRX) 513,517,517 
513 IF(DRX+ROE) 595,595,520 
517 IF(DRX-ROE) 520,595,595 
520 DRY=RY( J )-RY( I ) 

IF(DRY) 523,527,527 
523 IF(DRY+ROE) 595,595,530 
527 IF(DRY-ROE) 530,595,595 
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530 DRZ=RZ( J)-RZ(I ) 

IF(DRZ) 533,537,537 
533 IF(DRZ+RQE) 595,595,540 
537 IF(ORZ-ROE) 540,595,595 
5 40 DIST = DRX'^DRX+DRY*DRY+DRZ*DRZ 
IF(DIST-RaE2) 550,595,595 
550 DIST=SQRT(DIST) 

560 POT=EXP { PEXA+PEXB*DIST )-PPTC 
580 TPOT=TPQT+POT 

PPE( n=PPE( I )+BSAVE*PGT 
PPE{ J)=PPE( J)+TSAVE«POT 
595 CONTINUE 
600 CONTINUE 

200 DO 300 I=IQ,LL 

IF(LCUT( I ) ) 300, 205,300 
205 IP=I+1 

DO 295 J=IP,LL 
IF(LCUT(J)) 295,210,295 
210 DRX=RX( J )-RX( I ) 

IF(DRX) 213,217,217 
213 IF(DRX+RQEC) 295,295,220 
217 IF(DRX-ROEC) 220,295,295 
220 DRY=RY{ J )-RY{ I ) 

IF(DRY) 223,227,227 
223 IF{DRY+ROEC) 295,295,230 
227 IF(DRY-ROEC) 230,295,295 
230 DRZ = RZ( J)-RZ( I ) 

IF(DRZ) 233,237,237 
233 IF(DRZ+ROEC) 295,295,240 
237 IF(ORZ-RCEC) 240,295,295 
240 DIST=DRX^ORX+DRY^'=DRY+ORZY^DRZ 
IF(DIST-R0EC2) 250,295,295 
250 DIST=SQRT(DIST) 

IF(DJST-ROEA) 260,255,255 
255 IF{ niST-ROFB) 265,270,270 
260 POT=EXP( PV/i + EXR»::nT ST) 

GO TO 280 

265 POT = OIST*( DIST«{DI ST^CP3 + CP2 )+CPl )+CPO 
GO TO 280 

270 POT = EXP{ CGD1+CGB1=<^DIST ) -EXP ( CGD2+CGB2*D I ST ) 
280 TPOT=TPOT+POT 
SAVE=0. 5^P0T 
PPE( I )=PPE( I ) +SAVE 
PPE(J)=PPE( J)+SAVE 
295 CONTINUE 
300 CONTINUE 
RETURN 
END 



SUBROUTINE LOCAL 
C 

THIS SUBROUTINE CALCULATES THE TOTAL POTENTIAL ENERGY IN A 
SMALL VOLUME AROUND A VACANCY OR INTERSTITIAL. 

C 

COMMON/ CCM1/RX{ lOOC ) ,RY(1C00 ) ,RZ( 1000 ) , LCUTI 1000) , 
1LL,LD, ITYPE,NVAC 

CCMHON/COM7/PPTC,TPQT, PPE{ 1000 ) , TLP E , ROEL , ROEL 2 , NEW 
DIMENSION LCCAT{50C) 

K=1 

TLPE=0.0 

IFINEW.EQ. 1 )G0 TO 305 
8 GO TO { 10,20,20) ,ITYPE 
10 I=NVAC 
GO TO 200 
20 1 = 1 

200 DO 300 J=1,LD 

1F{ I.EQ. J) GO TO 250 
210 DRX=RX(J )-RX( I ) 

IF(DRX) 213,217,217 
213 IF(DRX+RGEL) 295,295,220 
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217 IF(DRX-ROEL) 220,295,295 
220 DRY=RY( J )-RY( I ) 

IF(DRY) 223,227,227 
223 IF(DRY+ROEL) 295,295,230 
227 IF(DRY-ROEL) 230,295,295 
230 DRZ=RZ( J )-RZ( I ) 

IF(DRZ) 233,237,237 
233 1F(DRZ+RQEL) 295,295,240 
237 IF(DRZ-ROEL) 240,295,295 
240 OIST=DRX-ORX+DRY*DRY+DRZ*ORZ 
ROEL2 = ROEL=^'ROEL 
IF( DIST-R0EL2) 250,295,295 
250 LOCAT(K)=J 
K = K + 1 

295 CONTINUE 
300 CONTINUE 
KF=K-1 

305 DO 310 K=1,KF 

TLPE=TLPE+PPE( LOCAT (K) ) 

310 CONTINUE 
NEW=1 
RETURN 
END 



SUBROUTINE PRINT 



THIS SUBROUTINE PRINTS THE HEADING OF ALL PERTINENT INFORMA- 
TION AT THE TOP OF EACH TIMESTEP PRINTOUT. 

C 

C0MM0N/C0M1/RX( 1000 ) , RY( 1000 ) , RZ( 1000 ) , LCUT( 1000) , 
1LL,LD, ITYPE,NVAC 

C0MM0N/C0M2/IHK 20 ) , I H2 ( 8 ) , I H5 (10 ) , I HB( 6 ) , I HT ( 6 ) , 
ITARGET (4) , TMAS , BULLET( 4 ) , BMAS, PL AN E , T EMP , T HERN 

CCMMCN/CCM3/RXI ( 1000) ,RYI ( lOOC ) ,RZI ( lOOC ) ,CVR,EVR, 

] NT.TTMF.DT.DTI , TLAY 

rnMMQ[\i/rnMA/rx - rv. i 7.<:ry.<;rv.sr7.ir)FEP-rii y.niv.ni? 

CCMM0N/C0M5/R0E , ROE 2 ) ROE M j E X A , E XB , PEX A i P E X B , FX A^ P FX A , 

1 IQ,TSAVE, BSAVE 

COMMON/ CCM8/R0E A , ROEB, ROEC , R0EC2 , CPO , CP 1 , CP2 , CP3 , 
1CFC,CF 1 , CF2,CGD1 ,CGD2, CGBl ,CGB2 ,CGF1 ,CGF2 

9710 FORMAT (40X, 10A4 , / , 2 8 X , 2C A4 , / ) 

9720 F0RHAT(9H TARGET - , 4A4 , ICHPR I M ARY - , 4A4 , 1 X , 1 4HL ATT I CE 
1 UNIT =,F7.4,4H ANG ) 

9730 F0RMAT(4X,6HMASS = , F7 . 2 , 1 3 X , 6HM AS S = , F7 . 2 , 9X , 14HL ATT I C 
IE TEMP =F5.2,7H DEG K,,18H THERMAL CUTOFF =,F5.2,3H F 
IV/ ) 

9740 F0RMAT(2H (,A4,8H) PLANE, ,18H 
1 F5.2,21HKEV, CRYSTAL SIZE ( 

1 ),, 4X, 16HVACANCY IN SITE , 

9741 F0RMAT(2H (,A4,8H) PLANE,, 18H 
1 F5.2,21HKEV, CRYSTAL SIZE ( 

1 ),, 4X, 15H1NTERSTITI AL ( - , F 5 . 2 , 2H , - , F 5. 2 , 2H , + , F 5. 2 , 

112H) FROM SITE , 14/) 

9742 FORMAT (2H (,A4,8H) PLANE,, 18H PRIMARY ENERGY =, 

1 F5.2,21HKEV, CRYSTAL SIZE ( ,I2,3H X , 1 2 , 3H X ,I2,3H 

1 ),, 4X, 20HRE PLACEMENT IN SITE , 14/) 

9750 F0RMAT(30H PRIMARY START POINT (LU) X =,F5.2,5H, Y =, 
1F5.2,5H, Z =,F5.2, 5X, 13, • LAYERS ARE FREE TO MOVE', 

1 10X,4HIQ =, 12/ ) 

9760 F0RMAT(12H POTENTIAL , 6A4, 3X , 5HP E XA= , F 9. 5 , 2 X , 5HP E XB= , 
1F9. 5,2X, 5HPFXA=, F9. 5) 

9765 FORMAT ( 12X, 6A4, 3X, 5HEXA = , F9. 5 , 2X , 5HEXB = , F9. 5 , 2 X , 5HF X 
lA =, F9. 5/ ) 

9770 FORMAT!' WHEN',F8.4,' < R <',F8.4,' THE MATCHING POTEN 

ITIAL PARAMETERS ARE',//,' CPO =',F10.3,', CPI =' 

1F10.3,', CP2 =',F10.3,', CP3 =',F10.3,/,' CFO =' 

1E10.3,', CFl =',E10,3,', CF2 =',E10.3,//) 

9780 FORMAT!' CUT-OFF AT',F5.2,', vn'HEN R > ',F6.3,' LU, MQR 
ISE POTENTIAL PARAMETERS A^E', 8A4,//,10X,' CGMl =', 



PRIMARY ENERGY =, 

, 12, 3H X , I2,3H X , 12, 3H 
14/ ) 

PRIMARY ENERGY =, 

, I2,3H X , I2,3H X , 12 ,3H 



lF8c4,', CGD2 =' ,E8o4, ' , CGBl =','^8.4, 
CGFl =',E8.4,', CGF2 =',F8.4,//) 



CGB2 =' 5 F8,^. , 
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9790 FORMATdOH TIMESTEP , 1 4 , 22 X , 6HDT I = , F5.3 t 5H LU, , 

122H ELAPSED TIME (SEC) =t E10.4,21H, LAST TIMESTEP WA 
IS =, E10.4/ ) 

WRITE ( 6,9710 IHS,IH1 
WRITE ( 6,9720) T ARGE T , BULL E T , C VR 
WRITE ( 6,9730) TMAS , BMAS , TEMP , THERM 
GO TO (401,402,403), I TYPE 

401 WRITE ( 6,9740) PL ANE , EVR , I X , I Y, I Z , NVAC 
GO TO 405 

402 WRITE ( 6,9741) P L ANE , E VR , I X , I Y , I Z , 0 1 X , 01 Y , D1 Z , N V AC 
GO TO 405 

403 WRITE ( 6,9742) PLANE , EVR , I X , I Y , I Z , NVAC 

405 WRITE ( 6,9750) R X I ( 1 ) , R Y 1 ( 1 ) , RZ I ( 1 ) , I LA Y , I Q 
WRITE ( 6,9760) I H B, P EXA , PEXB , PFXA 
WRITE ( 6,9765) I HT , EXA , EXB , FXA 

WRITE ( 6,9770) ROEA, ROE B , CPO , CPI , CP2 , CP3 , CFO , CF 1 , CF2 
WRITE ( 6,9780) ROEC , ROEB , I H2 , CGDl , CGD2 , C GB 1 , CGB2 , 
1CGF1,CGF2 

WRITE ( 6,9790) NT , DT I , T I HE , DT 
RETURN 
END 
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